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Abstract 


The first generation of stars, often called Population III (or Pop III), 
form from metal-free primordial gas at redshifts z ~ 30 and below. 


They dominate the cosmic star formation history until z ~ 15 — 20, at 


which point the formation of metal-enriched Pop II stars takes over. 


We review current theoretical models for the formation, properties and 


impact of Pop III stars, and discuss existing and future observational 


constraints. Key takeaways from this review include the following: 


Primordial gas is highly susceptible to fragmentation and Pop III 
stars form as members of small clusters with a logarithmically flat 
mass function. 

Feedback from massive Pop III stars plays a central role in regulating 
subsequent star formation, but major uncertainties remain regarding 
its immediate impact. 

In extreme conditions, supermassive Pop III stars can form, reaching 
masses of several 10° Mc. Their remnants may be the seeds of the 
supermassive black holes observed in high-redshift quasars. 

Direct observations of Pop III stars in the early Universe remain 
extremely challenging. Indirect constraints from the global 21 cm 
signal or gravitational waves are more promising. 

Stellar archeological surveys allow us to constrain both the low-mass 
and the high-mass ends of the Pop III mass distribution. Observations 
suggest that most massive Pop III stars end their lives as core-collapse 
supernovae rather than as pair-instability supernovae. 
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1. Introduction 


From studying the cosmic microwave background (CMB), we know our Universe started 
out very simple. It was by and large homogeneous and isotropic, with small fluctuations 
that can be described by linear perturbation theory. In stark contrast, the Universe today 
is highly structured on a vast range of length and mass scales. In the evolution towards 
increasing complexity, the formation of the first stars marks a primary transition phase of 
cosmic evolution. Their light ends the so-called ‘dark ages’, and they play a key role in 
cosmic metal enrichment and reionization, thereby shaping the Universe at large and the 
present-day galaxy population. The study of stellar birth in the early Universe is therefore 
important for many areas of modern astronomy and astrophysics but is also a relatively 
young field of research. Only with the advent of advanced numerical methods and powerful 
supercomputers did a comprehensive modeling of early star formation become feasible. As 
a consequence, there is still considerable debate about the physical processes that govern 
stellar birth at high redshifts and the overall properties of the first stars. This review aims 
at providing an overview of the current state of the field. 

The first generation of stars, the so-called Population III (or Pop III) build up from truly 
metal-free primordial gas. Second generation stars, sometimes termed early Pop II stars, 
form from material that has been enriched from the debris of the first stars. Unlike the very 
first stars, which have not yet been directly detected, members of the second generation 


2 Klessen & Glover 


have been found and characterized in surveys looking for extremely metal-poor stars in 
our Milky Way and neighboring satellite galaxies. In concert with observational data at 
high redshift, this allows us to constrain the properties of genuine Pop III stars. Initially, 
primordial star formation was proposed to be very simple, governed by well defined initial 
conditions provided by Gaussian fluctuations of the cosmic density field. These were thought 
to result in the build-up of solitary high-mass stars. This simple picture has undergone 
dramatic revisions. We now understand that fragmentation is a wide-spread phenomenon 
during first star formation, and that Pop III stars form as members of multiple stellar 
systems with separations as small as the distance between the Earth and the Sun. Studies 
that include radiative feedback, magnetic fields, dark matter annihilation, as well as the 
primordial streaming velocities add to this complexity. All of these processes are relevant 
and need to be included in any realistic model. There is agreement now that primordial 
star formation is just as dynamic and complicated as stellar birth at the present day. 

We structure our review as follows: First we provide in Section 2 an overview of all 
relevant physical processes that contribute to first star formation. This includes an intro- 
duction of the cosmological model that constitutes the basis of our analysis, together with 
the equations that govern the evolution of the cosmic fluid. We end this part with an 
account of the instabilities that lead to gravitational collapse and subsequently to stellar 
birth, and we stress the importance of the thermodynamic properties of the star-forming 
gas. In Section 3 we discuss the critical mass for collapse in the primordial Universe and 
describe the evolution of the cosmic star-formation rate density. We then discuss in detail 
the standard Pop III formation pathway in isolated pristine halos. We also account for 
more complex scenarios in which further physical processes add to the complexity of the 
problem. Finally, we turn our attention to the most extreme physical conditions that result 
in the formation of supermassive stars and lead to the supermassive black holes we observe 
in quasars at high redshift. The impact of stellar feedback is the focus of Section 4. We first 
look at the consequences of the radiative and the mechanical energy and momentum input 
from young massive stars on their immediate birth environment and then consider their 
influence on neighboring halos and on the larger-scale intergalactic gas. We also discuss 
chemical feedback of Pop III stars, which drives the transition to Pop II star formation 
and governs the early metal enrichment of the Universe. A critical assessment of possible 
observational probes of primordial star formation is presented in Section 5, covering mea- 
surements at high redshift as well as multi-messenger data from the local Universe. Finally, 
we conclude and summarize in Section 6. 

For further reading on primordial star formation, we refer to the reviews of Bromm & 
Larson (2004), Yoshida et al. (2012), and Bromm (2013), and to the book by Loeb (2010). 
A good overview with a specific focus on astrochemistry is provided by Glover (2005, 2013), 
and an account of important numerical aspects of high-redshift star formation is given by 
Greif (2015). For further reading on the transition to the second generations of stars and 
build up of the first galaxies we recommend Bromm & Yoshida (2011) or the textbook 
by Loeb & Furlanetto (2013). In order to calibrate our understanding of primordial star 
formation with what we know about stellar birth at present days, we also refer to the 
reviews by Mac Low & Klessen (2004), McKee & Ostriker (2007), Zinnecker & Yorke (2007), 
Krumholz (2015), or Klessen & Glover (2016) on different aspects of this subject. 
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Pop III stars: 
Population III stars 
form from truly 
metal-free 
primordial gas in the 
high-redshift 
Universe. 


Pop II stars: 
Population II stars 
build up from 
metal-enriched 
material. Early Pop 
II stars constitute 
the second stellar 
generation in cosmic 
history. 


Pop | stars: 
Population I stars, 
such as the Sun, are 
relatively young and 
metal rich. They are 
found in the thin 
disk of the Milky 
Way and other spiral 
galaxies. 


2. Important physical concepts 


In this Section, we introduce the basic physical concepts needed to understand the formation 
of the first stars. We begin with the cosmological model that forms the foundation of all 
further considerations, introduce the main equations governing the evolution of the cosmic 
fluid, and study the criteria for the onset of gravitational instability in the early Universe 
and consequently for the formation of the first stars. 


2.1. Cosmological model 


Throughout this review, our approach is based on the standard ACDM model (e.g. Peebles 
1993; Bullock & Boylan-Kolchin 2017) in which the Universe consists of matter in form of 
cold dark matter and baryons, radiation, as well as dark energy. In this model, the Hubble 
parameter H(a) = à/a, which describes the expansion of the Universe, evolves as 

H’ (a) Qr Qm QA 


SS at a + Bar with 1=0+94+9%, 1. 


where a(z) = (1 + z) `! is the cosmic scale factor and z is the redshift. At the present day, 
i.e. for a = 1, the densities of radiation, matter and dark energy with respect to the critical 
density are Q, 2 107?, Qm = 0.31 and QA = 0.69 (Planck Collaboration et al. 2020b). The 
parameter w relates the pressure and energy of the dark energy component (p = wpe’); here, 
we assume that it is constant in time and equal to -1, as appropriate for the classical cosmo- 
logical constant. The current value of the Hubble parameter Ho is ~ 70kms~'Mpc™!; for 
historical reasons, this is usually expressed as Ho = 100^ kms ! Mpc `! with h zz 0.7. Note 
that there remains some debate regarding the precise value of Ho — (for a comprehensive 
overview, see Valentino et al. 2021). 


2.2. Evolution of the cosmic fluid 


In order to understand primordial star formation we need to identify and characterize the 
first regions in the Universe that decouple from the cosmic expansion and contract under 
their own gravitational attraction. Observations of the CMB show that the Universe started 
out extremely simple. The initial density distribution was isotropic and almost perfectly 
homogeneous with spatial fluctuations, 6(x) = p(x)/(p) — 1, of order 107” on large scales 
with respect to the average background density (p). In the dark matter, these fluctuations 
are present on all scales and grow due to their own self-gravity. In the linear regime, the 
evolution of ó is governed by the equation 
0?6 06 


This corresponds to the famous Jeans (1902) equation describing the evolution of self- 
gravitating isothermal gas spheres, written in the limit of zero pressure, and with an ad- 
ditional drag term due to the expansion of the Universe. In the ACDM model, the initial 
density fluctuations are well characterized as a Gaussian random field, and so their statis- 
tical properties are completely determined by the power spectrum of the field, which has 
been inferred with high precision from the observed temperature fluctuations in the CMB 
(Planck Collaboration et al. 2020a). Consequently, the physical properties of the Universe 
at a redshift z ~ 1000, and hence the initial conditions for cosmic structure formation, are 
extremely well constrained. 
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2.3. Gravitational instability 


Whereas (cold) dark matter can be considered a pressure-less zero-temperature fluid, this is 
not the case for the baryonic component of the Universe. On large scales, the gravity of the 
dark matter dominates, and the evolution of the baryons is very similar to that of the dark 
matter. On small scales, the effects of gas pressure become increasingly important. In the 
classical Jeans (1902) stability analysis, we can identify a critical mass scale Mj separating 
the gravity-dominated and pressure-dominated regimes. Perturbations with M > MJ un- 
dergo gravitational collapse, whereas those with M < My; do not. In the simple case of an 
isothermal gas sphere, we have 


45/2 1 3⁄2 1⁄2 3 2 n -1/2 T 3/2 
= — | = 7 ES M p (=) — $ ; 
6 (G) = ane lcm-? (i) I 


where n is the number density of H nuclei and cs is the sound speed, which are related to 
T and p by p = umun and T = umuc?/kg. Information about the chemical makeup 
of the gas is encoded in the weight factor u = (1+ 4x) with x being the ratio of He to 
H atoms by number. For primordial gas, x = 0.079 and so u = 1.32. In this convention, 
atomic hydrogen at one particle per cubic centimeter corresponds to n = 1 cm^?, whereas 
fully molecular hydrogen with one particle per cm? yields n = 2cm ?. Note that the 
chemical weight factor u could simply be replaced by the mean molecular weight if we take 
the total particle number density in Equation 3, rather than just the number density of 
hydrogen atoms. Finally, we note that a similar expression to Equation 3 can be derived by 
considering the growth of plane-wave baryonic perturbations in the linear regime (see e.g. 
Peacock 1999). More generally, an expression differing only by a small numerical factor can 
be derived by comparing the gravitational and sound-crossing timescales or the gravitational 
and thermal energy of the perturbation. 

The critical mass for collapse can also be derived by considering the two competing 
timescales in the problem: 


37 YA : R : ; 
Tg = (5) (free-fall time), Tsouna = eg (sound-crossing time) . 4. 


s 


The dynamical or free-fall timescale, zm expresses the characteristic duration of gravitational 
collapse in the absence of pressure, whereas the sound crossing timescale, Tsound = R/Cs 
denotes the time it takes to communicate pressure gradients across a fluctuation of size 
R. If TE € Tsouna the system is unstable against contraction, and if Tsouna < Tg then 
pressure gradients are able to provide stability against gravitational attraction. Note that 
in equilibrium, these two timescales are the same. Note also that there is a third relevant 
timescale to consider, which is the age of the Universe. It is reasonably well approximated 
as the inverse of the Hubble parameter, m1 = 1/ H(z), tracing the cosmic expansion history. 
Even if T < Tsouna, if both numbers are larger than Ty there is not enough time for collapse 
to progress to sufficiently large densities for star formation to set in. 

For a collapsing sphere, we can estimate the associated accretion rate from the Jeans 
mass and the free-fall time 7g = V/31/32Gp via: 


s My c2 
M=f ` = f 5. 
poets. 


where f can take on values between ~ 1 and several tens, depending on the initial density 
profile and the ratio between M and M;. Note, in this simple approximation, the accretion 
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Jeans instability: 
Instability that leads 
to the collapse of 
self-gravitating 
isothermal spheres if 
gravity dominates 
over gas pressure. 


Primordial gas: 

Big Bang 
nucleosynthesis 
produces mostly H 
and “He with mass 
fractions of 0.76 and 
0.24, respectively, 
and trace amounts 
of 2H, ?He, and "Li 
at the level of 1073 
to 10719, 


Equation of state: 
It relates 
thermodynamic 
state variables such 
as pressure, density, 
temperature or 
internal energy, and 
is needed to turn the 
equations of 
hydrodynamics into 
a closed and thus 
solvable system. 


rate only depends on the gas temperature. For further discussion, see Whitworth & Sum- 
mers (1985) and references therein. Additional physical processes can also be accounted 
for in our definition of the Jeans mass, most easily by defining and using an effective sound 
speed in place of cs. For example, when the gas is turbulent on scales much smaller than 
the dynamical scales of interest (Chandrasekhar 1951a,b; von Weizsücker 1951), we can 
simply add the velocity dispersion c to the sound speed. A similar role has been ascribed 
to magnetic fields B, in which case we include a contribution from the Alfvén velocity, 
va = B/(8np)'/? (e.g. Federrath & Klessen 2012), leading to 


1 1/2 
Cs eff = E +o? + 3j : 6. 


We revisit these aspects later, in Section 3.3. 


2.4. Impact of thermodynamics 


The above considerations demonstrate the importance of the thermodynamic response of 
the evolving system (see also Section 3.2.1). As an illustration, look at the simplified case 
of the gas following an effective polytropic equation of state, 


po p^" with index ye —1-- dlnT/dlnp, 7. 


where the yer is the result of the competition between various heating and cooling mecha- 
nisms (e.g. Omukai et al. 2005; Klessen & Glover 2016). From Equation 3 we see that 


M; < pz ett 3) : 8. 


If yer > 4/3 the Jeans mass increases during the contraction and eventually becomes 
comparable with the mass of the system (including both dark matter and baryons). In 
this case pressure forces will stop further collapse. For adiabatic gas with yer = 5/3, as 
appropriate for monoatomic gas without internal degrees of freedom, this will happen long 
before stellar densities are reached. In order for stars to form, the gas must be able to 
radiate energy away during its collapse, so that yer remains below the critical value. 

The importance of cooling can be assessed by comparing the cooling time 


T _ 1 Tot kT 
RET 4y—1A(T,n)' 


where ntot is the total number density of particles, with the free-fall timescale, Tg. If 
Tcoo < Tg then the gas can cool rapidly, and collapse proceeds roughly on a free-fall 
timescale. However, if mm < 7Tcoo then the gas quickly becomes pressure supported, and 
the contraction slows down and proceeds quasi-statically for a duration of order of Tcoo1. 
As before, if 7 is the shortest timescale, we can consider the system as being stable. We 
note that the situation is typically more complicated than this simple analysis suggests, as 
the thermodynamic response of the gas may vary with time as density, temperature, and 
chemical composition evolve, and it may depend on location, specifically on the proximity 
to sources of stellar feedback, as discussed in Section 4. 


2.5. Instability of rotationally-supported systems 


In rotationally-supported systems, such as protostellar accretion disks or spiral galaxies, the 
criterion for gravitational instability takes a slightly different form, owing to the stabilizing 
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effect of shear. For infinitely thin disks, this was investigated by Toomre (1964), who derived 
the following criterion for instability: 


KCs 
<1. 10. 
TGM ^ 9 


Q= 
Here X and & are the surface density and epicyclic frequency, respectively. For systems in 
Keplerian rotation, we have k = Q, where Q is the rotational frequency (Kratter & Lodato 
2016). This approach can be extended to thick disks with multiple components (Rafikov 
2001; Elmegreen 2002) by introducing correction factors of order unity to Equation 10. 
There are two main pathways towards disk fragmentation. First, in the absence of 
accretion onto the system, Gammie (2001) argue that even an initially stable disk will 
become unstable if the cooling timescale, Tcoo1, is shorter than the orbital time scale Torbit = 
1/9. Second, in the presence of accretion from the surrounding gas envelope, if the mass 
load onto the disk exceeds its capability to transport material inwards, X increases beyond 
the critical value and the disk becomes unstable. As we see later, this latter scenario is 
commonly encountered when studying Pop III accretion disks. 


3. Pop III star formation 


We begin our discussion of Pop III star formation with a critical review of their birth 
environment, then turn our attention to the most likely formation pathway, and speculate 
about alternative scenarios. In the most extreme cases this can lead to the formation of 
supermassive stars as possible progenitors of supermassive black holes. 


3.1. Critical mass for collapse and cosmic star-formation rate density 


The formation of the first stars in the Universe occurs in regions where the cosmic fluid fulfils 
two conditions. First, it needs to decouple from the global expansion and begin to contract 
due to the self-gravity of dark matter. We call a region where this happens a dark matter 
halo. Second, the gas within the dark matter halo needs to be able to cool and go into 
run-away collapse to dramatically increase the baryon to dark-matter ratio and eventually 
reach stellar densities. As we discuss in more detail later, halos with virial temperatures 
above Tyir ~ 8000 K can cool initially via Lyman-a emission from atomic hydrogen, whereas 
lower mass halos with lower Tu, depend on cooling from H>. For historical reasons, the 
latter type of halo is often referred to as a ‘minihalo’. In the current ACDM paradigm, 
gravitationally bound objects form in a hierarchical fashion with smaller objects forming 
first, implying that minihalos are the first sites in which Pop III stars can potentially form. 


3.1.1. Simple models for the critical mass. Calculations of the growth of density perturba- 
tions in an expanding Universe (e.g. Barkana & Loeb 2001) show that substantial baryonic 
overdensities can only develop in halos with masses above a critical mass given by 


-1/2 2\ —3/5 3/2 
EU 3 Umh? Oph 1+z 


This mass scale is very similar to the Jeans mass (Equation 3) of fluctuations with a density 
close to the mean density of the Universe and a temperature that declines adiabatically as 
the Universe expands, and holds for z S 100 when gas is no longer thermally coupled to 
the CMB by Compton scattering. However, the development of a baryonic overdensity is 
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Toomre instability: 
Instability occuring 
in differentially 
rotating self- 
gravitating thin 
disks if gravity 
dominates over 
pressure and 
rotational shear. 


LW and UV photons: 
Lyman- Werner 
photons in the 
energy range 
11.2eV € hv « 
13.6eV can 
photodissociate H> 
molecules. UV 
photons with 

hv 2 13.6eV can 
photoionize H. 


Streaming velocity: 
Second-order cosmic 
perturbation theory 
predict a relative 
motion between dark 
mass and baryons, 
which decreases 
linearly as the 
Universe expands. 


a necessary but not sufficient condition for the formation of Pop III stars. In addition, as 
we have argued above, the gas must also be able to cool. This requirement yields a much 
higher critical mass. For instance, using a very simple model for the behavior of gas in a 
high-redshift minihalo Glover (2013) derives a value 


Mas ri 14 x 10°Mo ( 29 EE m pu 1+2\ 3°” 12 
ud © \ 0.143 1.32 10 ? i 


where we encounter again the chemical weight factor u introduced in Equation 3. 


3.1.2. Models with more physics included. Clearly, Equations 11 and 12 are substantial 
simplifications. For example, the latter assumes that the chemical evolution of the gas takes 
place at constant density and temperature, and neither equation accounts for the fact that 
conditions in the Universe change during the time that it takes for a perturbation to grow 
into the non-linear regime. This can lead to significant deviations from the prescriptions 
above (see e.g. Gnedin & Hui 1998 or Naoz & Barkana 2007 for more complex models). 
To capture the non-linear evolution and to account for the complex interplay between the 
different physical processes involved requires us to resort to numerical simulations. Processes 
that must be accounted for include the build-up of a cosmic Lyman-Werner (LW) or X-ray 
radiation background as star formation sets in (e.g. Gnedin 2000, see also Sections 4.3.2 and 
4.3.3) and the existence of relative streaming velocities between dark matter and baryons 
(Tseliakhovich & Hirata 2010; Tseliakhovich et al. 2011, see also Section 3.3.1). 

The latter results from second-order cosmological perturbation theory and therefore is 
often neglected in a purely linear analysis. Prior to recombination, baryons are tightly cou- 
pled to photons resulting in a standing acoustic wave pattern (Sunyaev & Zel’dovich 1970) 
and consequently in oscillations between baryons and dark matter with relative velocities 
of about 20 km s" and coherence lengths of of a few Mpc in comoving units at z © 1000 
(Silk 1968). After recombination, baryons are no longer tied to photons, their sound speed 
drops to ~ 6kms~", and the velocity with respect to the dark matter component becomes 
supersonic with Mach numbers of M ~ 5 (Tseliakhovich & Hirata 2010). As the Universe 
expands, the relative streaming velocity decays linearly and reaches ~ 1kms~ at z # 30, 
which is comparable to the virial velocity of the first halos to cool and collapse (Fialkov 
et al. 2012). Simulations that include streaming velocities suggest that their presence re- 
duces the baryon overdensity in low-mass halos, delays the onset of cooling, and leads to a 
larger critical mass for collapse to set in (Greif et al. 2011b; Stacy et al. 2011; Maio et al. 
2011; Naoz et al. 2012, 2013; O’Leary & McQuinn 2012; Latif et al. 2014a; Schauer et al. 
2017b; Nakazato et al. 2022). They may also have substantial impact on the resulting 21 cm 
emission (Fialkov et al. 2012; McQuinn & O'Leary 2012; Visbal et al. 2012). Increasing the 
LW background intensity also increases Merit, as discussed in more detail in Section 4.3.2. 

Altogether, there remains significant uncertainty about the most appropriate value of 
Merit with different simulations reaching different conclusions, depending on the numerical 
resolution adopted and the number of physical processes included (see, e.g. Yoshida et al. 
2003; Latif & Schleicher 2019; Kulkarni et al. 2021, in addition to the above). We illustrate 
some of these variations in Figure 1, where we plot predictions for Merit from Schauer et al. 
(2021), Kulkarni et al. (2021) and from models OHS and F13 in Hartwig et al. (2022) 
and Chen et al. (2022), which themselves are derived from the prescriptions of O'Shea & 
Norman (2008), Hummel et al. (2012), and Stacy et al. (2011) and from Fialkov et al. (2012, 
2013), respectively. For further details, see Appendix A.2 of Hartwig et al. (2022). The 


Klessen & Glover 


Va = 0.8+0.5 & Joi =0 
108 —— Schauer+21 
—— OHS === atomic cool. halo 
Bee — Fialkov+13 —'- Equation (18) 
107 ^ 77—— Kulami-21  :-- Naoz&Barkana07 
e 
z 
£ 106 
= 
10? 
10? 
10 15 20 25 30 35 40 
redshift 


Figure 1 


Critical mass for fluctuations in the cosmic fluid for collapse to set in. The solid lines show results 
from models that take the streaming velocity into account, adopted from Schauer et al. (2021), 
Kulkarni et al. (2021), and also models OHS and F13 from Hartwig et al. (2022) and Chen et al. 
(2022), which are derived from O’Shea & Norman (2008), Hummel et al. (2012), and Stacy et al. 
(2011) and from Fialkov et al. (2012, 2013), respectively. The fit from Schauer et al. (2021) is 
based on simulation data at z < 25 and we indicate our extrapolation of it to higher redshifts with 
a dotted line. The lines are computed for the most probable value of the streaming velocity, 

Ust = 0.8 (in units of the root-mean-squared value of the streaming velocity; see Schauer et al. 
2019) and the shading indicates the effect of varying this value by 0.5. For further comparison, 
we show an estimate of Merit in the absence of streaming based on our Equation 12 (dash-dotted 
line). The lower envelope (dotted line) shows the minimum halo mass required for a baryonic 
overdensity to develop, as computed by Naoz & Barkana (2007; see their Equations 17 and 18). 
The upper envelope (dashed line) is given by the critical mass for atomic cooling halos (see e.g. 
Equation 1 of Hummel et al. 2012 or Equation 15 of Hartwig et al. 2022). 


curves shown are computed in the absence of LW background radiation (J21 = 0), and for 
a streaming velocity vst = 0.8 (in units of the root-mean-squared value), which is the most 
likely value to be encountered in the Universe (Schauer et al. 2019, 2021). For a discussions 
of the large-scale impact of LW feedback and the associated uncertainties, we again refer 
to Section 4.3.2. For completeness, we also plot the estimate for Merit with vst = 0 given 
by Equation 12, the minimum mass required for the development of a baryonic overdensity 
computed by Naoz & Barkana (2007), and the critical mass for an atomic cooling halo (taken 
from Hummel et al. 2012 and Hartwig et al. 2022). These systems have a virial temperature 
of about 10* K, and their thermodynamic properties are dominated by Lyman-a cooling 
rather than H> cooling (Sutherland & Dopita 1993). Halos of this mass and above are 
able to cool and collapse even in the presence of a very strong LW radiation background 
(Oh & Haiman 2002) or a high streaming velocity (Schauer et al. 2019), and so this curve 
constitutes the upper envelope of all models for Merit considered here. 


3.1.3. Cosmic star-formation rate density. In the ACDM model, Pop III star formation 
begins at a redshift z > 30 in rare high-sigma fluctuations. The rate increases as collapse 


becomes possible in more and more halos and reaches a peak at redshifts z ~ 15 — 20. 
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Figure 2 


Star formation rate density for Pop III and Pop II stars at z > 5 from different studies, including 
numerical simulations (Johnson et al. 2013; de Souza et al. 2014; Xu et al. 2016; Jaacks et al. 
2018; Sarmento et al. 2019; Liu & Bromm 2020b) and semi-analytical models (Trenti & Stiavelli 
2009; Mebane et al. 2018; Visbal et al. 2020; Hartwig et al. 2022; Magg et al. 2022a). For 
comparison, we also provide the observational constraints from Madau & Dickinson (2014), 
Behroozi & Silk (2015), Finkelstein (2016), and Behroozi et al. (2019) with the gray shaded area 
indicating the proposed range. 


Although the overall cosmic star formation rate continues to increase (Madau et al. 2014), 
the rate at which metal-free Population III stars form declines again, because regions in 
the Universe that have not been enriched by supernova ejecta from massive stars become 
increasingly rare. This is the transition to the birth of (slightly) metal-enriched Population 
II stars, which we discuss in Section 4.4. Different models for the star formation rate density 
(SFRD) as function of redshift are depicted in Figure 2. The scatter gives a feeling for the 
current uncertainties in determining the onset and evolution of high-redshift star formation 
in both numerical simulations (Johnson et al. 2013; de Souza et al. 2014; Xu et al. 2016; 
Jaacks et al. 2018; Sarmento et al. 2019) and semi-analytical models (Trenti & Stiavelli 
2009; Mebane et al. 2018; Visbal et al. 2020; Hartwig et al. 2022). For comparison, we 
also provide the observational constraints from Madau & Dickinson (2014), Behroozi & 
Silk (2015), Finkelstein (2016), and Behroozi et al. (2019) with the total range indicated in 
gray, although we caution that the James Webb Space Telescope (JWST) will soon improve 
the situation for z >> 5. 
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3.2. Standard Pop III formation pathway 


We refer to stellar birth in pristine halos containing zero-metallicity gas, which have not 
been affected by stellar feedback from neighboring halos, as the standard Population III 
formation pathway. We briefly mention alternative scenarios in Section 3.3. 


3.2.1. Initial collapse phase. As discussed above, the ability of the gas in a halo to collapse 
and form stars depends on its ability to cool. All main cooling processes in zero metallicity 
gas are related to hydrogen, either in atomic or molecular form. At high temperatures, 
collisions can populate the excited electronic states of H which then de-excite by emitting 
Lyman series photons. This process is often referred to simply as Lyman-o cooling and 
is most efficient around temperatures of ~ 104K. To reach lower temperatures, molecular 
hydrogen is needed. The lightness of the H» molecule and its lack of a dipole moment 
conspire to render it an ineffective coolant at very low temperatures: its lowest energy 
radiative transition is the J = 2 — 0 transition in its vibrational ground-state, with an 
energy that corresponds to a temperature of ~ 512 K. The high-velocity tail in the thermal 
Maxwell-Boltzmann velocity distribution allows the gas to cool below this value, but only 
down to about 200 K (see Greif 2015). This is often termed the Pop III.1 formation pathway 
(e.g. McKee & Tan 2008; Clark et al. 2011a). The temperature can drop even further, if 
cooling by deuterated hydrogen (HD) takes over (Nagakura & Omukai 2005). HD has a 
non-zero dipole moment and its lowest energy transition is between the J = 1 to J = 0 
rotational levels, corresponding to a temperature of ~ 128 K. This is sometimes called Pop 
II.2 formation pathway, and in practice becomes relevant only in regions with enhanced 
fractional ionization, for example in very massive or in externally irradiated halos (see 
Sections 3.3 and 4.2). 

For a full account of primordial chemistry and the associated cooling and heating pro- 
cesses, we refer to the reviews by Glover (2005, 2013), or Bovino & Galli (2019). Here 
we focus on the most essential concepts. At low densities, H> has two main formation 
pathways. First, it can form by a two-stage reaction pathway involving the H ion as an 
intermediate step (McDowell 1961; Peebles & Dicke 1968). The reactions are: 


H+e > H +y, 13. 
H +H > H+e . 14. 
Second, there is also a contribution from a similar reaction pathway involving Hj as an 
intermediary molecule (Saslaw & Zipoy 1967): 
H+H" > Hi +y, 15. 
H} +H > H +H*. 16. 
Both reaction pathways require the gas to be partially ionized, and the amount of H> that 
can form via these routes is limited by the slow formation rates of the H^ and HZ ions 
and the recombination of the gas. Typically, the final molecular fraction is around 107°. 
At high particle densities above ~ 10? cm~?, a third process becomes important. It is the 
three-body reaction (Palla et al. 1983): 
H+H+H > H+H. 17. 
3 


As a result all atomic hydrogen is converted into Ha once particle densities of n zz 10'! cm7 
are reached. However, as collapse proceeds and the temperature exceeds values of ~ 2000 K 
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at densities of above 101? cm”? the Hz molecules become collisionally dissociated and the 
gas eventually turns atomic again. 
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The left side shows the relation between temperature T and hydrogen number density n during 
the initial collapse phase (top) and the corresponding fraction of molecular hydrogen (bottom). 
The data are taken from a resimulation of a halo from Schauer et al. (2021) provided by L. Prole. 
The right side (top) displays the relation between log1oT' and logion, corresponding to the 
effective equation of state of the gas, for different metallicities, ranging from a purely primordial 
composition (dark blue, Z = 0) to the solar value (red, Z = Ze) in steps of ten. The data are 
computed with a one-zone astrochemical model following the collapse of an isolated system (for 
full details, see Omukai et al. 2005, 2010). The dashed lines indicate constant Jeans mass for 
atomic gas (Equation 3). The right side (bottom) illustrates the main heating and cooling 
processes relevant for primordial gas. Cooling rates associated with hydrogen are depicted with 
solid lines, other cooling processes with dashed line, and the two main heating mechanisms with 
dotted lines. The data for the top right plot are provided by K. Omukai and the figure is 
reproduced with permission of ApJ. 


With the appropriate chemical network and the corresponding heating/cooling functions 
included in cosmic structure formation calculations we can follow the initial collapse phase 
from the cosmic mean up to the formation of the first hydrostatic core in the halo center. 
An example of this behavior is provided at the left side of Figure 3. The data are taken from 
a resimulation of a halo studied by Schauer et al. (2021) at a time just before the formation 
of the first protostar. We show the effective equation of state, i.e. the relation between 
temperature T' and number density n of hydrogen atoms (top), as well as the corresponding 
fraction of molecular hydrogen (bottom). The labels in the logo T — log) plot indicate 
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key phases of the initial collapse: (A) As the gas begins to flow into the potential well of 
the dark matter halo, it is compressionally heated to the virial temperature of the system. 
At the same time the H> fraction increases from < 107° to ~ 1073, which is sufficient 
for the gas to go into a run-away cooling phase (B) and brings it down to the minimum 
temperature of ~ 200K (C). As more gas flows into the halo center, eventually the potential 
becomes dominated by the gas rather than by dark matter. The contraction proceeds and 
the heat provided by PdV work begins to dominate again over cooling. As a consequence, 
the gas temperature rises to about 1000K at n ~ 10?cm ? (D). At this stage, three- 
body H> formation (Equation 17) becomes important and the gas quickly becomes fully 
molecular (E). Below densities n ~ 10!? cm”? the main cooling process is the (forbidden) 
ro-vibrational line emission from H>. However, at larger n two additional cooling processes 
become important (F). The first is collision-induced emission (CIE), which occurs when 
two molecules come close to each other. Van der Waals forces can then induce a temporary 
dipole which allows for efficient dipole emission during the interaction time interval (Omukai 
& Nishi 1998; Ripamonti & Abel 2004). The second one is associated with the collisional 
dissociation of H>, which sets in at temperatures around 2000 K and quickly dominates the 
overall cooling behavior. As the density increases further (G), more and more H> molecules 
are destroyed. This also implies that more hydrogen atoms become available again for 
H> formation through the three-body process and the associated energy release starts to 
dominate the overall heating rate over PdV work. This continues until the molecular gas 
is largely depleted and the collapse becomes almost adiabatic at n e 10?? cm ?. 

For comparison, we provide at the right side of Figure 3 (top) the relation between 
logioT' and logion from a one-zone astrochemical model describing the time-evolution of 
the central region during the initial collapse (Omukai & Nishi 1998) It essentially depicts 
the effective equation of state. We see that the polytropic index yer varies considerably 
across the different regimes introduced above. Data are taken from Omukai et al. (2005, 
2010) and describe the thermodynamic response of gas across a wide range of metallicities. 
We start with a purely primordial composition (Z — 0), indicated by the dark blue line, 
and show the results for logio Z/Ze = -6, -5,... up to the solar value (Z = Ze) in red. For 
comparison, we also show lines of constant Jeans mass following Equation (3) for atomic 
gas. 

For a more quantitative assessment of the most relevant heating and cooling processes 
during the initial collapse phase of gas in primordial halos, we also plot the corresponding 
rates as a function of the hydrogen number density, using data taken from the one-zone 
model of Glover & Savin (2009). Here, solid lines depict cooling processes associated with 
molecular hydrogen: ro-vibrational line emission from H> and HD, collision-induced emis- 
sion (CIE) of H> at high densities and collisional dissociation of H>. The cooling processes 
associated with other species are less relevant. Using dashed lines, we list cooling from Hi, 
LiH, Compton scattering, and H. In the absence of stellar feedback only two heating pro- 
cesses are important (dotted lines): PdV work dominates at densities below n ~ 10? cm ?, 
and the latent heat released by H> formation is relevant at higher densities. 

We note that over the roughly 18 orders of magnitude in density covered in Figure 3, the 
temperature of the primordial gas only varies by a factor of 25 or so. Overall the gas roughly 
exhibits an effective index of yer 7 1.08, which is close to the isothermal value of yer = 1. 
'This has important consequences for the level of fragmentation in primordial gas and it is 
essential to understand the dynamical evolution of the accretion disk that inevitably builds 
up around the central object (see Section 3.2.2). We also note that Figure 3 only considers 


www.annualreviews.org e The First Stars 


13 


14 


MH-1 | 5kpe (comoving) / | MH-I | 10 pe MH-1 | 200 AU 


MH-2 | 200 AU 


MH-3 | 200 AU 


Mone [Mo] 


log ny [cm] log ny [cm] log ny [cm?] 


e A 0 1 2 3 4 123 45 67 8 10 11 12 13 14 15 


Figure 4 


Structural and dynamical properties of three star-forming halos. The left side shows a sequence of 
zoom-ins focusing on the density-squared weighted hydrogen number density projected along the 
line of sight. The right plot provides cumulative mass as function of radius (top left) and 
corresponding number density profile (bottom left) in the inner regions of the halos (plus two 
additional systems not shown here) shortly before the onset of star formation. Specific angular 
momentum (top right), gas temperature (middle right), and turbulent velocity dispersion (bottom 
right) are plotted as functions of enclosed mass. The figure is adapted from Greif et al. (2011a) 
and reproduced with permission of ApJ. 


compressional and chemical heating. The situation may change if radiative feedback from 
newly formed stars (Section 4.2) or the possible energy input from turbulent dissipation or 
from dark matter annihilation is considered (Section 3.3.2). 


3.2.2. Disk formation and fragmentation. The first 3D simulations able to follow the initial 
cooling and collapse of the gas in high redshift minihalos became available around the year 
2000 (e.g. Bromm et al. 1999, 2002; Nakamura & Umemura 2001; Abel et al. 2000, 2002; 
Yoshida et al. 2003). Although these numerical simulations were fully three-dimensional, 
the halos considered were relatively round and so assuming spherical symmetry was a very 
good approximation during the early stages of collapse (see also Yoshida et al. 2006, 2008). 
These early calculations typically stopped when the object in the center reached hydro- 
gen number densities of n ~ 1019 cm 3, above which the computational timestep became 
prohibitively small. At this time the hydrostatic object in the center had a mass of only 
~ 107? Mo. At such an early time, the material that made it to the center carried very 
little angular momentum, and this object was surrounded by only a small disk-like struc- 
ture which was more strongly supported by pressure than by rotation. The authors of these 


Klessen & Glover 


studies suggested that this should be true for the entire protostellar accretion history, and 
hence argued that all the inflowing mass would end up in one single high-mass star (see also 
Tan & McKee 2004). Clearly this supposition needed to be tested, in particular, because 
the statement that primordial stars only form in isolation is in tension with present-day 
star formation, where fragmentation is ubiquitous and massive stars are typically found in 
clusters and aggregates (Lada & Lada 2003). 

The left side of Figure 4, adopted from a high-resolution simulation published by Greif 
et al. (2011a) ten years later, gives a visual impression of three star-forming halos at various 
spatial scales, eventually zooming in on the immediate physical environment where Pop III 
stars build up. The gas virializes on a scale of 5 kpc (comoving), followed by the runaway 
collapse of gas in the central 10 pc, where it becomes self-gravitating and decouples from 
the dark matter. In the final stages of the collapse, a fully molecular core forms on scales 
of ~ 200 AU. The right side of Figure 4 provides key structural and dynamical parameters, 
including cumulative mass and the corresponding radial density profile to the left, as well 
as specific angular momentum, gas temperature, and turbulent velocity dispersion as a 
function of enclosed mass to the right. Although the morphological parameters of the halos 
are very similar on large scales, they differ on small scales. They also have quite distinct 
dynamical properties. The collapse induces turbulence in the infalling gas, which leads 
to stochastic variations in the structural appearance of the central regions, demonstrating 
that small differences in the initial fluctuation spectrum can amplify during the non-linear 
collapse phase and lead to very different evolutionary pathways. For example, HD cooling 
became important in only two of the five halos modeled by Greif et al. (2011a), leading to 
lower temperatures in these halos and a different fragmentation behavior. This suggests 
that, similar to present-day star formation (McKee & Ostriker 2007; Klessen & Glover 
2016), stellar birth in the early Universe is also subject to large statistical variations with 
the outcome sensitively depending on the details of the initial conditions and environmental 
parameters as well as on various complex non-linear feedback loops. 

Current simulations of Pop III star formation acknowledge this complexity and attempt 
to follow the formation and dynamical evolution of the accretion disk that builds up around 
the central object, and ideally cover the entire accretion history, until stellar feedback re- 
moves the remaining gas and the process of stellar birth is completed. These studies demon- 
strate that primordial accretion disks are highly prone to fragmentation. They suggest that 
the standard pathway of Pop III star formation leads to a stellar cluster with a wide range 
of masses rather than the build-up of one single high-mass object. This is illustrated in 
Figure 5, taken from Clark et al. (2011b), showing the early evolution of a typical Pop II 
accretion disk. It fragments and builds up a multiple system of four protostars within only 
about one hundred years after the formation of the first object. 

Disk fragmentation across a wide range of spatial and temporal scales is reported in 
essentially all current studies (e.g. by Machida et al. 2008a,b; Turk et al. 2009; Clark et al. 
2011b; Greif et al. 2011b, 2012; Smith et al. 2011, 2012a; Dopcke et al. 2013; Susa 2013; 
Susa et al. 2014; Vorobyov et al. 2013; Stacy & Bromm 2013; Stacy et al. 2016; Hirano 
et al. 2014; Hosokawa et al. 2016; Takahashi & Omukai 2017; Hirano & Bromm 2018; Susa 
2019; Wollenberg et al. 2020; Sugimura et al. 2020; Sharda et al. 2021; Jaura et al. 2022; 
Chiaki & Yoshida 2022; Prole et al. 2022b, which is by no means an exhaustive list). They 
predict the formation of a cluster of multiple Pop III (proto)stars, which grow in mass at 
rates of M z 107? Mo yr 1 with possible brief periods of accretion with rates as high as 


1 


a few times 102 Mo yr 1. The reason for the high susceptibility of primordial accretion 
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Figure 5 


Formation and evolution of the disk around a Pop III protostar, adopted from Clark et al. 
(2011b). The left side gives a visual impression of the build-up of a two-arm spiral structure 
caused by gravitational instability which subsequently fragments into four protostars within 

110 yr. The right side shows the radial surface density and temperature profile at different times. 
It also demonstrates that the disk material is almost fully molecular (ng, /n = 0.5) and Toomre 
unstable (Q « 1). The figure is reproduced with permission of Science. 


disks to fragmentation is always the same: For typical halo conditions in the primordial 
Universe the mass load onto the disk from the infalling envelope exceeds its capability to 
transport material inwards by gravitational or magnetoviscous torques. As a consequence 
massive spiral arms build up and speed up the inward transport. Often this is not enough, 
and the arms become non-linear and interact with each other, leading to run-away collapse 
in the interaction regions, as clearly visible in Figure 5. The disk is fragmenting and forms 
new protostars. As this happens, a smaller accretion disk forms around the new object, 
which by itself may become unstable and fragment to form additional protostars. For very 
quiescent initial conditions with low levels of turbulence, this secondary process may be the 
prevalent form of fragmentation (Susa 2019). As the global accretion disk gains more mass 
and grows in size by accretion of higher angular momentum material, the Toomre unstable 
region moves further out. Consequently, fragmentation and formation of new protostars 
occurs at larger and larger radii as the evolution progresses. 


3.2.3. Pop III IMF and multiplicity. As argued in the previous Section, it is almost impos- 
sible to avoid the fragmentation of primordial accretion disks. The key question then is, 
what happens to the fragments during their subsequent dynamical evolution? There are 
three possible outcomes: (1) A fragment grows in mass and survives to become a proper 
star (or substellar object), or (2) it dies by getting swallowed by the central object, or (3) it 
merges with another fragment, in which case the same question applies to the new object. 
A related question is, if the fragment survives, what sets its final stellar mass? Again, 
there are three answers possible: (1) The system simply runs out of mass, leaving behind 
relatively massive stars with interesting implications for their detectability at high redshifts 
(Section 5.1). (2) Stellar feedback removes material from the system and terminates any 
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subsequent accretion. This is the most complex case and can lead to a wide range of differ- 
ent results, with the currently available simulations being highly inconclusive, as discussed 
in Section 4.2.1. (3) The fragment gets ejected from the disk by dynamical interactions and 
its growth is stopped prematurely. This typically leads to low-mass objects, some of which 
could potentially have survived until the present day and might be detectable in stellar 
archeological surveys (see the discussion in Section 5.2). 

Let us consider these processes in more detail. As material flows through the disk 
towards the center, it first encounters the Hill volume of protostars that lie further out, 
and so it preferentially gets swallowed by these objects rather than by those close to the 
center. This reduction of the mass growth rate is more pronounced at smaller radii and can 
lead to the complete starvation of the primary object in the very center (Peters et al. 2010; 
Girichidis et al. 2011). In simple binary systems, this process results in the two constituent 
stars having roughly equal masses (see also Kratter & Matzner 2006). In the more complex 
situation of Pop HI star formation, with ejections, mergers, and absorption by the central 
object, all of which are unpredictable and highly stochastic processes, we expect a wide 
spectrum of stellar masses. 

When considering the question of what fraction of fragments merge or become swallowed 
by the central object and how many survive, various studies (e.g. Greif et al. 2012; Smith 
et al. 2012a; Stacy & Bromm 2013) indicate that roughly 2/3 of the fragments quickly 
disappear again, and about 1/3 remain. These numbers should be taken with caution, 
as none of the simulations covers the entire duration of disk evolution, with the highest 
resolution models being able to cover the least amount of time. It therefore could be that 
protostars that are counted as ejected will eventually fall back again and become accreted 
then. Also those remaining in the disk may still merge at a later stage. In addition, the 
result also depends on how mergers are numerically implemented (Wollenberg et al. 2020). 
Still, as long as fragmentation in the disk continues and the appearance and disappearance 
of new protostars is an ongoing process, adopting a ratio of roughly 1:2 of survivors and 
mergers is a good estimate. It has also been suggested that the number Nkrag of fragments, 
survivors as well as mergers, increases with time t as Nfrag cx ips (Susa 2019). 

Even if Pop III stars form in multiple systems across a wide range of masses, we can 
nevertheless look for the most bound pairs and investigate their properties. Doing so, Stacy 
& Bromm (2013) and Stacy et al. (2016) report a binary fraction of ~ 35%, with semi-major 
axes as large as 3000 AU and with a wide distribution of orbital periods, ranging from ~ 30 
to 30000 years, with the smallest period being determined by the minimum numerical 
resolution achieved. They also find that the distribution of mass ratios is relatively flat (as 
also suggested by Susa 2019). With these data, it is possible to build large ensembles of 
Pop III binary systems for the long-term integration in dedicated N-body simulations or 
semi-analytic models (e.g. Liu et al. 2021; Santoliquido et al. 2021). 

Concerning the resulting distribution of stellar masses, the initial mass function (IMF) 
of Pop III stars, the ubiquity of fragmentation and the stochasticity of the processes involved 
lead to a wide range of stellar masses, with all current models suggesting values from the 
substellar regime up to several hundred solar masses. As mass is the most important 
parameter in stellar evolution theory (Kippenhahn et al. 2012), predicting the Pop III 
IMF is a central aspect of many studies. In Figure 6 we present the results of some of 
these simulations to give an indication of the diversity of the current state-of-the-art. We 
select three high-resolution models without stellar feedback: (1) model psink = 1075 gcm”? 
from Prole et al. (2022b); (2) model 005001 from Wollenberg et al. (2020); and (3) the 
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Figure 6 


Normalized differential (top) and cumulative (bottom) stellar mass distribution obtained from six 
different studies of Pop III star formation. To represent the current state-of-the art, we select 
three high-resolution models without stellar feedback (Prole et al. 2022a: model with 

Psink = 10-9 gcm-?, Wollenberg et al. 2020: model a005001, Stacy & Bromm 2013: all models) 
and three simulations with radiative feedback included (Jaura et al. 2022: model T.RTP, Hirano 
et al. 2015: model IIL.1, Hirano et al. 2014: all models). Despite the large variation in range, the 
predicted IMFs are rather flat compared to the present-day IMF. 


combined model from Stacy & Bromm (2013). Models (1) and (2) reach very high spatial 
resolution, but only cover a small fraction of the total accretion time in the halo. Model 
(3) is not so well resolved, but covers more time. The effect of resolution is noticeable in 
the resulting mass spectrum. We also chose three simulations with radiative feedback from 
newly formed stars: (4) T-RTP from Jaura et al. (2022); (5) IIL.1 from Hirano et al. (2015); 
and (6) the combined models from Hirano et al. (2014). Model (4) has similar resolution 
and time coverage as the first two examples (1) and (2) without feedback, and it predicts 
a similar mass spectrum. Models (5) and (6) have somewhat lower resolution, but cover 
considerably more time until the protostellar accretion has stopped. Whereas (1) — (4) are 
fully three-dimensional simulations, (5) and (6) have been performed in two dimensions 
assuming axisymmetry. These calculations predict considerably larger masses. Altogether, 
the reported stellar mass spectra vary enormously, which is likely a consequence in part of 
the difference resolutions and periods simulated in the different calculations. However, it 
is a common feature that they all are approximately logarithmically flat. T'his is in stark 
contrast to the present-day IMF, which shows a clear peak around 0.2 — 0.3 Mo followed 
by a steep power-law fall-off (Kroupa 2002; Chabrier 2003). For this reason, the IMF of 
primordial stars is often called top-heavy in comparison. 


Klessen & Glover 


We conclude that our understanding of the IMF of primordial stars is still quite limited. 
The existing models all have specific shortcomings, with none being able to reach sufficiently 
high resolution, cover enough time, and properly include all relevant physical processes. As 
indicated in Figure 6 we find that models that include stellar feedback on average predict 
larger stellar masses than those without. In addition the outcome strongly depends on the 
spatial resolution achieved, on the time span covered, and on the details of the numerical 
implementation. High-resolution simulations tend to yield smaller stellar masses, because 
they are better able to resolve disk fragmentation (see, e.g., the discussion in Stacy et al. 
2016 or Prole et al. 2022b), but typically only cover a small fraction of the full accretion 
history of the halo. This raises the question of whether the remaining gas reservoir is 
mostly used to form new (low-mass) objects, or whether it is largely consumed by existing 
protostars, which then continue to grow in mass. Similarly, simulations covering a large 
fraction of the halo collapse timescale do not resolve the inner accretion disk well, and 
therefore exhibit lower levels of fragmentation. Hence, they are biased towards higher-mass 
stars. Furthermore, two-dimensional simulations tend to report less fragmentation than full 
three-dimensional ones. Fragmentation is also influenced by additional physical processes, 
which we discuss in Section 3.3 below. 


3.2.4. Intermediate conclusions. Altogether, we can draw the following conclusions from 
the discussion so far: First, once a halo decouples from the cosmic expansion, the question 
of whether or not it forms stars depends strongly on its ability to cool. Second, as gas flows 
in and assembles in an accretion disk it is highly susceptible to fragmentation. Consequently, 
the formation of binaries or higher-order multiple stellar systems is the norm of Pop III star 
formation rather than the exception. Third, when focusing on binary stars and the most 
bound pairs in a hierarchical system, we expect a wide range of separations, a flat mass ratio 
distribution, and a roughly thermal spread of orbital eccentricities. Fourth, the resulting 
stellar mass spectrum is very wide, ranging from low-mass stars, which potentially could 
have survived until the present day with interesting implications for stellar archaeological 
surveys (Section 5.2.1), to very massive stars, which produce copious amounts of UV and 
LW photons in the early Universe. Overall the resulting mass spectrum is likely to be rather 
flat (in the logarithm of mass) and therefore top-heavy compared to the present-day values. 
Fifth, developing a better understanding the impact of stellar feedback is one of the key 
challenges of current research into Pop III star formation. Specifically, radiative feedback 
is very likely to reduce the number of stars within a star-forming halo and lead to larger 
stellar masses (see also Section 4.2). However, it can also potentially increase the level of 
(large-scale) fragmentation in externally irradiated halos (Section 4.3). 


3.3. Alternative Pop III star-formation pathways 


Here we touch upon alternative Pop III star formation scenarios and speculate about the 
impact of physical processes that are often neglected in numerical and analytical models. 
We group our discussion into processes that can promote fragmentation in the primordial 
gas and those that can reduce the level of fragmentation. We defer the discussion of the 
most extreme physical conditions, which might lead to the formation of supermassive stars 
and the seeds of supermassive black holes in the Universe, to Section 3.4, and we note that 
a detailed analysis of the impact of stellar feedback is the focus of Section 4. We also note in 
passing that adopting different cosmological models compared to standard ACDM (Section 


www.annualreviews.org e The First Stars 


19 


20 


2.1) can further modify the picture. It is well understood, for example, that in a warm dark 
matter scenario star formation happens later and in larger halos (Maio & Viel 2015; Dayal 
et al. 2015; Magg et al. 2016; Mocz et al. 2020). A similar result has been reported for fuzzy 
dark matter models (e.g. Mocz et al. 2019). However, it remains to be seen whether these 
global changes actually influence the properties of the individual stars that form, and we 
do not follow up further on this topic here. 


3.3.1. Enhancing fragmentation. Several numerical studies indicate that fragmentation also 
occurs on larger scales in the halo, on scales of the star-forming cloud as a whole (Turk et al. 
2009; Stacy et al. 2010; Clark et al. 2011b) rather than just at the level of the central disk, 
as discussed in Section 3.2.2. It is typically associated with higher levels of turbulence in 
the halo gas. The initial turbulent flows that are always present with subsonic or transsonic 
velocities, get amplified during gravitational collapse and induce density fluctuations that 
can go into run-away growth in their own right. This is very similar to the turbulence-driven 
mode of star formation that is dominant at the present day (Scalo & Elmegreen 2004; 
Elmegreen & Scalo 2004; Mac Low & Klessen 2004; McKee & Ostriker 2007; Krumholz 
2015; Klessen & Glover 2016). There are two systems, for which this effect is thought to 
be particularly important: atomic cooling halos and halos that are subject to large relative 
streaming velocities between baryons and dark matter. 


External irradiation and atomic cooling halos. Atomic cooling halos are halos with 
virial temperatures high enough to allow cooling by atomic hydrogen, i.e. Tyir > 8000 K. In 
general, we expect Pop III stars to form in a halo once its virial temperature exceeds a few 
thousand K (see Section 3.1) and hence that most atomic cooling halos will be associated 
with metal-enriched gas. However, in the presence of a sufficiently strong LW radiation field, 
H2 cooling and star formation can be suppressed in halos with Tyir < 8000 K (see e.g. Oh & 
Haiman 2002; Agarwal et al. 2019, and also Section 4.3.2). In that case, cooling will only 
get underway once the halo has grown to the point that Tyir > 8000 K. Then the gas will 
start to cool quasi-isothermally via Lyman-a emission. During this initial period of cooling, 
the high gas temperature keeps the critical mass for run-away collapse high (see Section 
3.1). This implies that once collapse becomes possible, the associated infall velocities and 
rates, which we can infer from Equation 5, are larger than the ones discussed in Section 
3.2.1. The high gas temperature also keeps the gas partially ionized, enabling it to form H» 
rapidly once it becomes dense enough to shield itself from the external LW radiation field, 
resulting in a transition to efficient H> cooling and a rapid temperature drop (Oh & Haiman 
2002). The combination of large velocities and low temperatures means that the turbulence 
associated with the inflow motion is likely to become trans- and supersonic (Greif et al. 
2008; Wise & Abel 2007a; Wise et al. 2008). The flow is characterized by cold streams that 
bring dense material rapidly to the center, where it can efficiently fragment and form stars. 
It has been suggested that this so-called Pop II.2 mode of primordial star formation leads 
to a different IMF (McKee & Tan 2008; Clark et al. 2011b; Maio et al. 2011; Stacy et al. 
2011), but overall the results are not fully conclusive. 


Streaming velocity between baryons and dark matter. High levels of turbulence are 
also expected in regions of large relative streaming velocity between baryons and dark 
matter (Section 3.1.2). Simulations that include this effect suggest that it reduces the gas 
overdensity in low-mass halos, delays the onset of cooling, and leads to a larger critical 
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mass for collapse to set in (see e.g. Greif et al. 2011b; Stacy et al. 2011; Maio et al. 2011; 
Naoz et al. 2012, 2013; O’Leary & McQuinn 2012; Latif et al. 2014a; Schauer et al. 2017b, 
2020, 2021). It may also have substantial impact on the resulting 21cm emission (Fialkov 
et al. 2012; McQuinn & O’Leary 2012; Visbal et al. 2012, see also Section 5.1.3). Gas in 
halos that are subject to large streaming velocities is also more turbulent than gas in more 
quiescent systems, and so we expect more fragmentation and a bias towards smaller stellar 
masses (Clark et al. 2008). However, this process has not yet been modeled with sufficient 
resolution, and so no reliable predictions of its impact on the IMF exist. 


3.3.2. Reducing fragmentation. Besides stellar feedback, which we address in detail in Sec- 
tion 4, two additional physical processes have been suggested to reduce the level of frag- 
mentation and thus to influence stellar IMF and multiplicity. These are magnetic fields and 
dark matter annihilation, both of which we discuss below. 


Magnetic fields. The presence of dynamically important magnetic fields could significantly 
alter the picture presented so far. We know that the current Universe is highly magnetized 
on all scales (Beck et al. 1996) and that this influences the birth of stars and the evolution 
of the interstellar medium.’ The properties of the magnetic fields observed today are well 
explained by a combination of small-scale and large-scale dynamo processes (Brandenburg 
& Subramanian 2005). In contrast, our knowledge of magnetic fields at high redshifts is 
very sparse. Theoretical models predict that magnetic fields could be produced in vari- 
ous ways, for example via the Biermann battery (Biermann 1950), the Weibel instability 
(Lazar et al. 2009; Medvedev et al. 2004), or thermal plasma fluctuations (Schlickeiser & 
Shukla 2003). Other theories place their origin in cosmological phase transitions or during 
inflation (Sigl et al. 1997; Grasso & Rubinstein 2001; Banerjee & Jedamzik 2003; Widrow 
et al. 2012). The resulting fields are thought to be orders of magnitudes too weak to have 
any dynamical impact, and so magnetohydrodynamic effects have often been neglected in 
numerical simulations of primordial star formation (however, see the analytic models of 
Pudritz & Silk 1989; Tan & McKee 2004; Silk & Langer 2006). 

This situation has changed with the realization that the small-scale turbulent dynamo 
can efficiently amplify even extremely small primordial seed fields to the saturation level 
(Kulsrud et al. 1997), and that this process is very fast, acting on timescales much shorter 
than the free-fall time. An analytic treatment is possible in terms of the Kazantsev model 
(Kazantsev 1968; Subramanian 1998; Schober et al. 2012a,b). This describes how the twist- 
ing, stretching, and folding of field lines in turbulent magnetized flows leads to exponential 
growth of the field. The amplification timescale is comparable to the eddy-turnover time 
on the viscous or resistive length scale, depending on the value of the magnetic Prandtl 
number (the ratio of the kinematic viscosity to the magnetic diffusivity). Once backreac- 
tions become important, the growth rate slows down, and saturation is reached within a 
few large-scale eddy-turnover times (Schekochihin et al. 2004; Schober et al. 2015; Liu et al. 
2022). Depending on the properties of the turbulent flow the magnetic energy density at 
saturation is thought to lie between 1% and a few 10% of the kinetic energy density with a 
field topology that is highly tangled (Federrath et al. 2011; Seta & Federrath 2021). 

Magnetic fields with that strength can strongly affect the evolution of protostellar ac- 
cretion disks. If the field has a strong polodial component, it can efficiently remove angular 


lFor the extreme viewpoint of magnetically mediated star formation, see Shu et al. (1987). 
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momentum from the star-forming gas and reduce its level of fragmentation (Machida et al. 
2008a,c,b; Machida & Doi 2013; Bovino et al. 2013; Latif et al. 2013, 2014b; Hirano & 
Machida 2022; Saad et al. 2022), which influences the resulting IMF (Turk et al. 2011; 
Peters et al. 2014; Sharda et al. 2021). A highly ordered small-scale field can also drive pro- 
tostellar jets and outflows (Machida et al. 2006; Sadanari et al. 2021). Finally, the presence 
of a dynamically significant magnetic field can also change the rotational properties of Pop 
III stars (Machida et al. 2007; Stacy et al. 2013), which in turn has consequences for their 
expected lifetimes and overall luminosity (see Section 4.2). One important caveat, however, 
is that many studies of the impact of the magnetic field assume that the field is highly- 
ordered on small scales. This is not true initially for a field amplified by the small-scale 
turbulent dynamo, and although we expect the field to become more ordered over time, 
it remains unclear how rapidly this occurs, with recent simulations yielding contradictory 
results (Sharda et al. 2021; Prole et al. 2022a; Stacy et al. 2022; Saad et al. 2022) 

Altogether, we expect Pop. III clusters to have fewer members with somewhat higher 
masses than predicted by purely hydrodynamic simulations, such as discussed in Section 3.2. 
However, the details depend very much on the adopted field topology, with initially smooth 
and homogeneous fields leading to more magnetic braking then highly tangled turbulent 
fields (e.g. Seifried et al. 2012; Kuruwita & Federrath 2019) as expected to result from the 
small-scale turbulent dynamo at work (e.g. Seta & Federrath 2020). And so the questions 
of how magnetic fields influence the overall star-formation process in primordial gas and 
how they affect the resulting IMF remain subject to very active research. 


Dark matter annihilation. Despite its importance for cosmic evolution and structure 
formation, the true physical nature of dark matter is still unknown. Many models intro- 
duce a new class of weakly interacting massive particles (WIMPs), as they naturally occur 
in supersymmetry theories (e.g. Jungman et al. 1996). The lightest supersymmetric par- 
ticle is expected to be stable and to have properties consistent with the phenomenological 
requirements on dark matter (Bertone et al. 2005). If these particles are self-annihilating, 
they will act as an additional source of heating. 

In most environments the dark matter density is too low for this to be significant. 
However, this may be different in the very centers of star-forming halos in the early Universe. 
Here, the collapse of the baryons may lead to adiabatic contraction of the dark matter halo 
(Blumenthal et al. 1986), increasing its central density by several orders of magnitude. As 
the annihilation rate scales quadratically with density, the corresponding energy input and 
ionization rate may become large enough to influence gas dynamics. Spolyar et al. (2008) 
and Freese et al. (2009) suggest that this process may overcome the cooling provided by 
H>, and speculate that this could halt gravitational collapse and lead to the formation of 
so-called dark stars. If dark matter particles also scatter weakly on baryons, these dark 
stars would be stable for a long time without ever becoming dense or hot enough to initiate 
nuclear fusion. They would be much larger and more massive than normal Pop. III stars, 
with sizes of a few AU, lower surface temperatures, and higher luminosities (Freese et al. 
2008; Iocco 2008; Iocco et al. 2008; Yoon et al. 2008; Hirano et al. 2011). 

There are several problems with this scenario. First, it is not clear whether collapse 
stalls once the energy input from dark matter becomes comparable to the cooling rate. Ri- 
pamonti et al. (2010) argue that this is not the case because the larger heating rate catalyzes 
further formation of H2 and is compensated by the corresponding larger cooling rate. In ad- 
dition, chemical cooling due to H> collisional dissociation can help to balance the additional 
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heating after only a moderate increase in the gas temperature (Smith et al. 2012a). Second, 
the implicit assumption of perfect alignment between dark matter cusp and gas collapse is 
most likely violated in realistic star formation conditions. Three-dimensional simulations 
(Stacy et al. 2012, 2014) clearly demonstrate that the presence of non-axisymmetric pertur- 
bations leads to a separation between dark matter cusp and collapsing gas, rendering the 
annihilation energy input insignificant for dark star formation. However, it is still possible 
that dark matter annihilation influences the dynamics of the accretion disk in the halo 
center and that the energy input associated with this process leads to a suppression of disk 
fragmentation (Smith et al. 2012a). We conclude that dark matter annihilation may be able 
to reduce the multiplicity of metal-free stars and increase their overall mass, but we also 
note that the existing studies are still premature and it is too early for a reliable assessment 
of the impact of this process on the IMF of Pop III stars. 


3.4. Supermassive stars and black holes 


The existence of extremely bright quasars at redshifts z > 7 (see e.g. Fan et al. 2006b; Wu 
et al. 2015; Banados et al. 2018; Wang et al. 2021) implies the presence of supermassive 
black holes (SMBH) with masses of M, ~ 10? Mc and above. This finding is in tension 
with physical models in which the seeds of SMBH start out light and then grow gradually 
by Eddington-limited accretion. The standard Pop III star-formation pathway discussed in 
Section 3.2 produces black holes at the end stage of stellar evolution with M, < 102 Mo 
(see also Section 4.3.4). These small seeds do not have enough time to grow to the observed 
large masses over the age of the Universe (~ 0.7 Gyr at z = 7, assuming standard ACDM 
cosmology). In addition, the assumption of persistent Eddington-limited accretion is itself 
questionable. Stellar feedback from Pop III stars is very efficient at removing gas from their 
birth sites, leaving little to be accreted by their black hole remnants, and the motion of 
these black holes within their parent minihalos can also be substantial, further reducing 
their accretion rate (Smith et al. 2018). Therefore, more extreme formation scenarios need 
to be considered. We briefly review a few of the most important models here and refer the 
reader to Woods et al. (2019), Latif & Schleicher (2019), or Inayoshi et al. (2020) for more 
comprehensive reviews. 


Need for high initial infall rates. If we assume spherical symmetry and if we furthermore 
assume that the opacity of the material falling onto the black hole is given by Thompson 
scattering on free electrons, then the released radiation must not exceed the Eddington 
luminosity, 


Lgaa = 4nMepmpc/or ; 18. 


where mp is the proton mass, p is the mean molecular weight and reflects the chemical 
composition of the gas (see also the discussion below Equation 3), c is the speed of light, 
and or is the Thompson scattering cross section. Otherwise radiation pressure would be 
too strong to allow for accretion. This can be converted into a characteristic timescale 
by comparing with the available energy reservoir in the system coming from the rest mass 
energy of the black hole, rgaa = Mec? /Leaa = ovc/ (4v ump), which is independent of the 
black hole mass. For typical conditions in primordial gas, "pa ~ 0.4 Gyr. The time taken 
for a black hole to grow via Eddington-limited accretion from an initial mass M., to a final 
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mass M, r is given by 


TEdd € M. d 
T. = In i I 19. 
feaa 1 — € E 


where e is the accretion efficiency, typically ~ 0.1 for accretion through a thin disk (Shakura 
& Sunyaev 1973), and where fraa expresses the fraction of time for which the black hole is 
able to accrete at the full Eddington rate (for a comprehensive review, see Inayoshi et al. 
2020). To build a SMBH of 10? Mo at z = 7 starting from an initial seed mass of 10? Mo 
requires T. 7 0.8 Gyr, which is longer than the age of the Universe at this redshift, even if 
we make the unrealistic assumption that feaa = 1 throughout. Reducing feaa makes this 
discrepancy even larger. In principle, accretion at a rate faster than the Eddington limit 
is possible if the accretion flow strongly deviates from spherical symmetry and material 
gets delivered to the black hole in a highly filamentary fashion, or alternatively, if the 
accretion disk is radiatively inefficient or the accreting envelope emits anisotropically (e.g. 
Mayer & Bonoli 2019). However, the inferred small quasar duty cycles would require rates 
considerably above the Eddington limit during the periods of accretion. This seems unlikely, 
and so the preferred scenario for the formation of SMBH at high redshifts is to start out 
with more massive seeds, i.e. to increase Me, to a much larger value. 

To build such a massive seed object requires very high accretion rates, which in turn 
requires extreme environmental conditions. Clearly, the halo in which this happens needs 
to be massive enough to contain a sufficient amount of gas. To form a seed with a mass of 
~ 10° Mo, we therefore need a halo with a mass M >> 10’ Mo, unless we assume that an 
improbably large fraction of the gas ends up in the seed black hole. This is much larger than 
the critical halo mass required for Pop III star formation and hence points to scenarios in 
which collapse and star formation is delayed in some fashion. Inayoshi et al. (2020) review 
the different models that have been suggested to explain this, ranging from irradiation of 
the halo by a high flux of LW photons (e.g. Agarwal et al. 2012) to dynamical heating of 
the gas by repeated major mergers (e.g. Mayer & Bonoli 2019). 

Once we have physical conditions that allow a large amount of gas to flow into the 


center of an appropriately massive halo at rates of M ~ 0.1 Mọ yr 1 


or higher, possibly 
reaching up to 1000 Mo yr `! (Zwick et al. 2023), the next question is whether this gas feeds 
the growth of a single object, or whether the gas fragments and forms multiple objects. We 


briefly consider both possibilities below. 


Supermassive stars. lf we assume that only a single object forms — a scenario often 
referred to as direct collapse — then we can ask about the properties of this object and 
about the maximum permitted mass. These questions can be addressed with detailed 
stellar evolution calculations that include a treatment of accretion (e.g. Stahler et al. 1986; 
Behrend & Maeder 2001; Omukai & Palla 2003; Hosokawa & Omukai 2009; Haemmerlé 
et al. 2016). They show that stable supermassive stars (SMS) are possible and can reach 
maximum masses of up to several 10° Mg (Umeda et al. 2016; Woods et al. 2017; Haemmerlé 
et al. 2018), but always stay below 10° Mo for the accretion rates reachable in atomic cooling 
halos (Haemmerlé et al. 2019; Haemmerlé 2020, 2021a,b). This value is determined by the 
onset of the general-relativistic (GR) instability, which drives the object into collapse to 
form a black hole with the same mass (see Appenzeller & Fricke 1972a,b or Fuller et al. 
1986 for early models, or Haemmerlé et al. 2021 for a more recent account). In Figure 7 we 
report estimates of the maximum mass from different sets of stellar structure and evolution 
calculations for different accretion rates and from simplified analytic models. The stellar 
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Figure 7 


Maximum mass reachable by supermassive stars accreting at different rates until the onset of the 
GR instability triggers collapse to a black hole of the same mass. The colored symbols depict the 
final masses in the studies of Umeda et al. (2016), Woods et al. (2017), and Haemmerlé et al. 
(2018) for accretion rates up to M — 10 Mo yr-1. We also report results from Haemmerlé et al. 
(2019) for 100 Mọ yr”! and 1000 Mç yr-!, but note that the calculations terminated before the 
GR instability sets in, and so the expected final masses are higher than shown here, as indicated 
by the tentative extension of the accretion path (vertical dotted lines). We also indicate the region 
forbidden by the GR instability as derived from analytic estimates (see Haemmerlé 2020, 2021a). 
The tilted lines show the evolution in M vs. M space for the free-fall collapse of gas spheres with 
different initial density profiles, characterized by the contrast p/p between the density of the 
infalling gas and the mean density of the accreting SMS, following the mass-radius relation for 
high accretion rates (Equation 23). The figure is inspired by Haemmerlé et al. (2021). 


structure calculations indicate that primordial SMS evolve as red supergiant protostars 
(Hosokawa et al. 2012a, 2013; Haemmerlé et al. 2018), with extended radii that follow the 
relation given by Equation 23 and surface temperatures of only ~ 5000 K. Their internal 
structure consists of a convective core, a radiative zone containing most of the stellar mass, 
and a convective envelope that covers a dominant fraction of the photospheric radius. To a 
good approximation, these structures can be described as hylotropes (Begelman et al. 2008; 
Begelman 2010), in particular for accretion rates > 10 Mọ yr !. Although very luminous, 
their low surface temperatures do not allow them to emit large amounts of ionizing photons. 
Consequently, they are not able to create extended HII regions (Section 4.2.1), which might 
limit the overall mass growth or affect star formation in neighboring halos (Section 4.3.1). 


Dense clusters. As argued in Section 3.2.2, it is difficult to prevent primordial gas from 
fragmenting. This is well established in the standard Pop III formation scenario, but also 
appears to hold in highly irradiated atomic cooling halos where H> is suppressed (e.g. 
Agarwal et al. 2012; Sugimura et al. 2014; Agarwal & Khochfar 2015; Latif et al. 2015, 2020; 
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Regan & Downes 2018). It becomes even harder to prevent fragmentation in the presence 
of metals or dust, owing to the additional cooling channels that they provide (Omukai 
et al. 2008; Latif et al. 2016; Chon & Omukai 2020). Therefore, rather than forming a 
single star, it is plausible that the rapid accretion flows discussed here instead feed the 
growth of a dense cluster of objects, each accreting at its own pace. In this scenario, SMBH 
seeds result from runaway collisions between stars in this dense environment (e.g. Portegies 
Zwart & McMillan 2000, 2002; Sesana et al. 2005; Devecchi & Volonteri 2009; Devecchi 
et al. 2012; Katz et al. 2015; Sakurai et al. 2017; Reinoso et al. 2018, 2020; Escala 2021; 
Vergara et al. 2021) or between their black hole remnants (Davies et al. 2011; Lupi et al. 
2014; Antonini et al. 2019; Kroupa et al. 2020). Although many models have treated this 
as a purely stellar dynamical problem, the difficulty of disrupting these high accretion flows 
means that in practice one should simultaneously take gas dynamics and stellar dynamics 
into account. In this context, Davies et al. (2011) and Lupi et al. (2014) argue that the 
inflow of gas into an interacting cluster of stellar mass black holes is needed to steepen the 
gravitational potential and make mergers more likely than three-body ejections. Focusing 
on embedded star clusters, Boekholt et al. (2018) demonstrate that the combination of 
accretion and collisions can lead to masses up to ~ 10° Mọ. The semi-analytic models of 
Tagawa et al. (2020) consider the growth of a supermassive object via stellar bombardment 
in the presence of gas. Other models study the impact of different accretion prescriptions 
(Das et al. 2021b), or of mass loss occurring in mergers (Alister Seguel et al. 2020) or 
associated with stellar winds (Das et al. 2021a). The impact of varying the metallicity of 
the gas has also been investigated (Chon & Omukai 2020; Schleicher et al. 2022). 

The emerging picture is that fragmentation leads to the formation of a dense and deeply 
embedded cluster in which competitive accretion of individual cluster members in concert 
with frequent merger events results in the run-away growth of a small number of objects. 
Depending on the environmental conditions and on the detailed implementation of the 
physical processes considered, masses of order of M, ~ 10° Mo are easily within reach, 
and it is reasonable to speculate that the internal structure of these very massive run-away 
objects is similar to that of the SMS discussed above. This is supported by calculations 
adopting highly time-varying accretion rates (Woods et al. 2021a,b) and considering SMS 
mergers as particularly extreme forms of accretion spikes. They collapse into massive black 
holes once this run-away growth phase ends or once they reach the mass limit for the general 
relativistic instability. 


4. Feedback from Pop III and transition to Pop II 


The formation of Pop III stars is associated with a range of different feedback processes 
that affect the gas around them. Radiative and mechanical feedback in various forms 
influences the star formation efficiency of the gas on both small and large scales, and 
chemical feedback, i.e. the enrichment of the Universe with metals from the first supernovae, 
drives the transition from Pop III to Pop II star formation. In this section, we review the 
most important feedback processes. We begin with radiation (Section 4.1), and successively 
include other forms of stellar feedback, focussing first on their effects on gas close to the 
stars (Section 4.2) and later on their effects on much larger scales (Section 4.3). We also 
briefly review the physics of the Pop III to Pop II transition and the most important open 
questions associated with this (Section 4.4). 
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4.1. Radiation from Pop III stars 


Here we summarize the properties of the radiation produced by primordial stars, which sets 
the stage for the discussion of the impact of radiative feedback on the star-forming cloud 
itself and on neighboring halos. 


4.1.1. Accretion luminosity. The first form of feedback to become important in a star- 
forming minihalo is the accretion luminosity generated by gas accreting onto newly-formed 
protostars. In principle, this represents a considerable reservoir of energy, owing to the high 
accretion rates typical in systems forming Pop III stars. A protostar with a mass M. and 
radius R. that is accreting gas at a rate M will produce an accretion luminosity 


GM.M 
ge 


Lace = " 
R. 


20. 
where o is a dimensionless efficiency factor that depends on the geometry of the accretion 
flow (Stacy et al. 2016), with high resolution simulations (e.g. Wollenberg et al. 2020; Jaura 
et al. 2022) suggesting values of a close to unity. Accretion rates onto Pop III protostars can 
vary substantially from star to star, but at early times, values as high as 10 ?—10 2 Mo yr "1 
are commonly encountered (Section 3.2.2). For a low mass pre-main sequence Pop III star 
with M, = 1Mọ and R, = 20Ro (Omukai & Palla 2003), this corresponds to accretion 
luminosities in the range Lace ~ 2000 — 20000 Lo, which is much higher than the main 
sequence luminosities of these same stars. Accretion onto more massive Pop III protostars 
produces even higher luminosities, with Lace scaling approximately as Lace cx MLT for 
fixed M (Smith et al. 2012a). 

Despite this, the effect of this radiation on the surrounding gas is relatively modest. 
Rapidly accreting Pop III protostars have photospheric temperatures of at most ~ 6000 K 
(Stahler et al. 1986b; Omukai & Palla 2003) and hence radiate most of their energy at visible 
and infrared wavelengths where the continuum opacity of metal-free gas is very small (Mayer 
& Duschl 2005). Therefore, only a tiny fraction of the radiated energy is absorbed by the 
gas in the minihalo, with most escaping into the intergalactic medium (IGM). The heating 
that this provides to the minihalo gas moderately reduces its propensity to fragment (Smith 
et al. 2011) but otherwise has little impact on its evolution. 


4.1.2. Stellar luminosity. At later times, once the newly formed Pop III stars have joined 
the main sequence (MS), their intrinsic luminosity becomes more important than their ac- 
cretion luminosity, specifically when talking about high-mass stars above ~ 10 Mo. Even 
on the MS, if the stars continue to accrete the geometry of the flow has strong impact on 
the internal structure (Hosokawa et al. 2012b). There are two extreme cases that we can 
consider. Accretion through a geometrically thin disk-like structure delivers fresh mate- 
rial with relatively low entropy. It settles onto the stellar surface with the same specific 
entropy as the photosphere, which allows the Pop III star to retain a relatively compact 
configuration. The total stellar luminosity L. relates to the radius R. and effective surface 
temperature T. as 

L, = Anosp RT , 21. 


with esp being the Stefan-Boltzmann constant. Consequently, a star subjected to cold disk 
accretion is relatively hot and compact. If the accretion flow is spherically symmetric, the 
infalling material retains the entropy produced in the accretion shock front, and the star 
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Figure 8 


Summary of key properties of Pop III stars in the mass range M. = 1.7 Mo to 500 Mo. To the 
left we show the stellar loci in the Hertzsprung-Russell diagram (logio L vs. logio T) across their 
MS and post-MS phases. For comparison, we also show the ZAMS of solar metallicity stars in the 
range of 0.8 Mç to 120 Mọ. To the right, we present the numbers for stellar lifetime, stellar 
radius, total luminosity, and effective surface temperature, each as a function of mass. The symbol 
size indicates the evolutionary stage, with the largest one indicating the ZAMS and the smallest 
the post-MS giant phase, and with a step size that corresponds to 1096 of the total stellar lifetime. 
The solid blue line illustrates the mass-radius relation for massive main-sequence Pop III stars 
(Equation 22). The underlying data are taken from Murphy et al. (2021a,b) and Martinet et al. 
(in prep.) for zero metallicity, and from (from Ekstróm et al. 2012) for solar metallicity. More 
details are provided in the Appendix. 


responds to the delivery of high-entropy material by increasing its radius. This implies that 
the surface temperature must drop, and so stars experiencing spherical hot accretion tend 
to be relatively cool and bloated (Hosokawa et al. 2012b). 

The exact values of the stellar radius R», effective surface temperature T, and total 
luminosity L = Lacc-- L, vary with time and strongly depend on the instantaneous accretion 
rate M and mass M. (Hosokawa & Omukai 2009; Hosokawa et al. 2010, 2012a; Haemmerlé 
et al. 2018). However, once the stars have finished accreting, they settle onto a well-defined 
main sequence. Pop III stars on this main sequence are more compact than their present-day 
counterparts, with a mass-radius relation at the high mass end given roughly by 


0.6 
R. ~ 03 Ro [V] : 22. 


They are therefore somewhat hotter than their present-day counterparts (Maeder & Meynet 
2012). We illustrate key properties of Pop III stars for masses ranging from M. = 1.7 Mo 
to 500M in Figure 8. The values are based on the stellar evolution models computed by 
Murphy et al. (2021a,b) and Martinet et al. (in prep.) using the Geneva code (Eggenberger 
et al. 2008). To the left we plot the location of the stars in the logio L-logio T diagram 
during the main-sequence (MS) and post-MS evolution. For comparison, we also provide 


the zero-age main sequence (ZAMS, see e.g. Kippenhahn et al. 2012) of solar metallicity stars 
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(from Ekström et al. 2012) in the mass range from 0.8 Mọ to 120 Mọ. This demonstrates 
that primordial and present-day stars of the same initial mass have similar total luminosity 
during the MS evolution, before mass loss becomes important for solar metallicity stars. 
However, due to their much higher compactness, Ts is considerably larger. To the right 
we visualize the stellar lifetime, as well as the stellar radius (with the blue line justifying 
Equation 22), the total luminosity, and the effective surface temperature as function of mass 
at various evolutionary phases. The largest symbol depicts the stars at the ZAMS and the 
smallest one in their giant phase. Altogether we provide these properties in steps of 10% 
of the total stellar lifetime. We note that different stellar evolutionary models and codes 
(e.g. SEVN, see Spera et al. 2022) lead to very similar numbers. Further information can 
be found in the Appendix together with tables of the fluxes of ionizing and non-ionizing 
photons, and the total number of photons produced per baryon, including also population- 
averaged values. 

In summary, the existing models strongly suggest that high-mass Pop III stars are very 
compact with high photospheric temperatures (Omukai & Palla 2003), and therefore radiate 
a large fraction of their energy at ultraviolet wavelengths. Once Pop III stars with masses 
greater than about 10 — 20 Mç have formed, they become powerful sources of ionizing 
photons. These photons can have a profound effect on the surrounding gas, as we explore 
below in Sections 4.2 and 4.3. 


4.1.3. Influence of rapid accretion on the stellar spectrum. The simple picture sketched 
above assumes that once a Pop III star reaches the main sequence, it remains there in a 
compact state. This is a reasonable assumption if the accretion rate onto the star is less 
than ~ 10 ? Mo yr ! (e.g. Hosokawa & Omukai 2009; Hosokawa et al. 2010; Maeder & 
Meynet 2012). However, if the accretion rate onto the star becomes larger, this triggers the 
expansion of the outer envelope of the star to a red supergiant state, in which the stellar 
radius is given by (Hosokawa et al. 2012a, 2013; Haemmerlé et al. 2018) 


23. 


M. ) 


EN 2 


This can be analytically derived by taking the Eddington luminosity, which is realistic 
for photon-dominated massive stars (e.g. Kippenhahn et al. 2012), adopting an effective 
surface temperature of T. ~ 5000 K, which is appropriate for red supergiants and which 
turns out to be roughly constant across a wide range of masses, and inserting both values 
into Equation 21. For comparison, the radius of a zero-age main sequence Pop III star with 
M. = 100 Mo is R. ~ 5 Re according to Equation 22. The intrinsic luminosity of the star 
is not strongly affected by this change in its structure, and Equation 21 implies that the 
increase of R, leads to a substantial drop of Ts down to values of ~ 5000 K, remaining 
relatively constant across a wide range of masses. This leads to a massive drop in the 
ionizing photon output of the star (Hosokawa et al. 2016). If the Pop III star is located 


l or above onto the star can 


in an environment where an accretion rate of 107? Mo yr” 
be maintained for an extended period, the star never develops an HII region. Therefore, 
photoionization feedback never has the opportunity to interfere with accretion onto the 
star, and it can grow to become extremely massive (as discussed in Section 3.4). 
Furthermore, even if the mean accretion rate onto the star is less than 107? Mo yr "1, 
it can still grow to reach a supergiant state if the accretion is highly episodic, with periods 


where the accretion rate exceeds the critical value (Vorobyov et al. 2013). In this case, a 
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HII region: 

Bubbles of ionized 
gas surrounding 
high-mass stars or 
massive star clusters 
caused by radiative 
stellar feedback. 
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period of rapid accretion triggers the expansion of the outer envelope of the star, but once 
the accretion rate falls below the critical value, the star contracts back to its original size. 
This process can happen repeatedly during the growth of the star (Hosokawa et al. 2016) 
and leads to the repeated appearance and disappearance of the surrounding HII region, in 
a fashion similar to the “flickering” found in simulations (e.g. Peters et al. 2010; Galván- 
Madrid et al. 2011) and observations (e.g. De Pree et al. 2015) of present-day HII regions. 

How common this behavior is in real systems depends on how often young Pop III 
stars encounter such high accretion rates. This is difficult to predict robustly, because it 
depends not only on the larger-scale environment in the halo, but also on the details of the 
disk fragmentation, the proportion of fragments that merge, etc. Simulations of this yield 
numerically converged results only when run with very high resolution (Prole et al. 2022b) 
and as yet very few such simulations have been run for long enough to also capture the 
formation of massive stars (see the discussion in Section 3.2.3). 


4.2. Small-scale impact 


In this section, we summarize the impact that feedback from Pop III stars has on gas close 
to the stars, which we here take to mean gas within the same dark matter halo. The impact 
of Pop III stellar feedback on much larger scales is discussed in Section 4.3 below. 


4.2.1. Photoionization. As argued above, high-mass Pop III stars are thought to produce 
copious amounts of ionizing radiation. Consequently, we begin our discussion of stellar 
feedback by investigating the impact of photons with energies above 13.6eV. 


Impact on surroundings. The ionizing photons produced by a massive Pop III star are 
readily absorbed in the surrounding gas, leading to the formation of an HII region. The 
evolution of this HII region once it breaks out of the dense protostellar accretion disk has 
been modelled by a number of different groups (e.g. Whalen et al. 2004; Wise & Abel 
2008; Stacy et al. 2016; Hosokawa et al. 2016) and is reasonably well understood. Because 
of the high densities found in the disk, the gas there resists photoionization and so the 
HII region initially expands in a bipolar fashion above and below the disk midplane, as 
illustrated in Figure 9. The gas within the HII region is hot — temperatures of 1— 2 x 10^ K 
are typical (Whalen et al. 2004) — and highly over-pressured compared to the surrounding 
neutral gas. It therefore expands outwards, driving a shock in front of it and leaving behind 
a low density cavity. The expansion of the HII region reduces the flow of gas onto the 
protostellar accretion disk, since gas can now reach the disk only by moving inwards close 
to the midplane. In addition, the disk itself begins to be photoevaporated by the radiation 
from the star (McKee & Tan 2008). Together, these effects lead to the loss of the disk 
and the termination of accretion onto the star after a few 10“ yr (Hosokawa et al. 2016; 
Sugimura et al. 2020). 

'The final fate of the HII region depends upon its location. In a minihalo, the fact that 
the sound speed of the photoionized gas exceeds the escape velocity of the minihalo means 
that the HII region will continue to expand for as long as the Pop III star continues to emit 
ionizing photons. For minihalos with masses close to Merit, the end result is that a large 
fraction of the gas mass is driven out beyond the virial radius (Kitayama et al. 2004; Whalen 
et al. 2004; Alvarez et al. 2006; Abel et al. 2007). In more massive minihalos, although the 
ionized gas is too warm to remain gravitationally bound to the minihalo, the larger size 
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t= 1.83e+4 yrs 
M. = 264.6 Ms 


temperature [K] density [/cc] velocity [km/s] 
( left-half ) ( right-half ) (lower values) 
3.0e+4 1.0e+10 10.0 
5335 3.2e+08 5.62 
948.7 1.0e+07 3.16 
168.7 3.2e+05 1.78 
30.0 1.0e+04 1.0 


t=1.17e+4 
M. 


yrs 


velocity 
(higher values) 
M : 100 [km/s] 


Figure 9 


Illustration of the bipolar expansion of HII regions around massive Pop III stars. Five examples 
from the simulations of Hosokawa et al. (2016) are shown. In each case, the Pop III star is located 
at the center of the double panel and the temperature, density and velocity are shown in the plane 
containing the polar axis of the disk. Velocities below 10km s^! are denoted by colour-coded 
arrows of fixed size; those above 10km s^! are denoted by white arrows, with lengths varying in 
proportion to the velocity. All five HII regions are shown at the time at which their polar 
extension first reaches 10* AU. This time is indicated in the figure for each HII region, along with 
the mass of the Pop III star at this moment. Figure is from Hosokawa et al. (2016) and 
reproduced here with permission of the authors and ApJ. 


of the system means that a single HII region is unlikely to drive out all gas. Finally, in 
protogalaxies with Vese ~ 20km s! or more, the ionized gas remains gravitationally bound 
to the protogalaxy, which therefore retains most of the gas. For a similar reason, the time- 
averaged escape fraction of ionizing photons from the HII region, i.e. the fraction of photons 
that can can travel into the IGM, also depends strongly on the location of the HII region. In 
small minihalos, it can reach values of 5096 or more (see e.g. Kitayama et al. 2004; Alvarez 
et al. 2006; Wise et al. 2014; Schauer et al. 2015), while in more massive galaxies, a broader 
range of values is possible (see e.g. Paardekooper et al. 2013; Dayal & Ferrara 2018). 


HII region trapping. An additional complication has been highlighted by Jaura et al. 
(2022). Most numerical studies of the growth of Pop III HII regions implicitly assume that 
the HII region has already escaped from the dense protostellar accretion disk by choosing 
an injection region for the ionizing photons that is larger than the disk scale height near the 
star.” Jaura et al. (2022) show that if one instead injects the photons at the location of the 


?The studies by Stacy et al. (2016) and Sugimura et al. (2020) do not make this assumption. 
Instead, they each adopt a different sub-grid model to treat the behavior of the HII region on 
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star itself, then the HII region remains trapped in the dense gas and does not break out of 
the disk. The physical reason for this is that at the high densities found close to the star, 
the Strómgren radius is much smaller than the Bondi-Parker radius at which the thermal 
pressure of the ionized gas first starts to dominate over the gravitational attraction of the 
star. Therefore, the ionized gas remains gravitationally bound to the star, preventing it from 
flowing outwards and escaping (Keto 2007). This result holds regardless of whether or not 
one accounts for the radiation pressure exerted on the gas by the ionizing photons, as this is 
too weak to overcome gravity close to the star. However, Jaura et al. (2022) also note that 
there are two important physical processes that are not included in their model that might 
change this picture. First, their simulations do not include a magnetic field and hence are 
unable to capture the formation of a magnetically powered jet or outflow. If such an outflow 
forms, then this will remove dense gas from the vicinity of the star, potentially allowing the 
HII region to break out from the disk (McKee & Tan 2008). Second, the high optical depth 
of the accretion disk to Lyman-a radiation means that any of the photons produced by the 
star or in the HII region will scatter repeatedly before escaping. The resulting radiation 
pressure is orders of magnitude larger than the direct radiation pressure of the ionizing 
photons and hence may be sufficient to overcome the gravitational attraction of the star 
(McKee & Tan 2008; Jaura et al. 2022). Modelling this will require a careful treatment of 
the interaction between the radiation field and the gas in the dense disk environment, and 
should be a priority for future simulations. 


4.2.2. Photodissociation. As well as producing a large number of ionizing photons, massive 
Pop III stars also emit many photons in the Lyman-Werner bands of H> (11.2 — 13.6eV). 
Because of the crucial importance of H> cooling in Pop III star formation, photodissociation 
of H» by these photons is potentially an important form of negative feedback. 

At high densities (e.g. in the immediate vicinity of the disk), the photodissociation 
region (PDR) created by these photons is restricted to a narrow region just outside of the 
HII region. This is because as dissociating photons penetrate into the gas ahead of the 
ionization front (I-front), they build up a dense layer of atomic hydrogen. Once the column 


density of this layer exceeds Ny ~ 1074 cm”? 


, absorption of LW photons in the damping 
wings of the atomic hydrogen Lyman series lines starts to be highly effective (Wolcott- 
Green & Haiman 2011). This shields the gas ahead of the dense atomic layer and limits 
its maximum column density to ~ 1025 cm ?. In the density regime where three-body H> 
formation is important (n » 10! cm ?; see Figure 3), the maximum spatial extent of this 
layer is only Lmax ~ 1014 ni cm z 10 ni AU, where nii is the gas density in units of 
10!! cem-?. Photodissociation of Hz in this narrow PDR leads to a dramatic increase in 
its temperature, driven primarily by H> formation heating (Susa 2013). Unfortunately, the 
impact of this temperature increase remains unclear. Some authors (Susa 2013; Susa et al. 
2014; Stacy et al. 2016) report that it drives a pressure wave outwards that reduces the 
gas density near the star and dramatically reduces its accretion rate even in the absence 
of photoionization feedback. On the other hand, Hosokawa et al. (2016) find it to have a 
much more limited effect and argue that photoionization is the dominant process. However, 
none of these studies accounts for the shielding of H2 by atomic hydrogen and hence they 


all overestimate the amount of dense gas affected by photodissociation feedback. 


scales that are not resolved in their calculation. However, neither of these subgrid models properly 
captures the behavior of the HII region near the star, as discussed in Jaura et al. (2022). 
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At lower densities, shielding of H> by atomic hydrogen becomes ineffective, as the 
required shell thickness becomes implausibly large. For example, at a density of ng = 
10° cm~?, producing a column density of 107° cm”? requires a shell thickness of ~ 7000 AU, 
several times larger than the entire physical extent of the region with mean density 
> 105 em" 
to the HII region and instead can expand to large distances within ~ 10^ yr. Gas in the 
PDR is not strongly heated by photodissociation, which is therefore ineffective at driving 


(see Figure 4). Therefore, the PDR is no longer restricted to the gas close 


gas out of the minihalo. However, the loss of the main coolant from the gas leads to it 
heating up to a temperature of 5000 — 10000 K as it collapses (Susa et al. 2014). Over time, 
this heating reduces the flow of gas to the center of the minihalo, limiting the amount of 
mass available for accretion even in the absence of effective photoionization feedback. 


4.2.3. Outflows. Protostellar outflows, both in the form of highly-collimated jets and also 
wide-angle winds, are a ubiquitous feature of present-day star-forming regions (Bally 2016). 
These outflows are magnetocentrifugal in nature. They are powered by the rotation of the 
protostellar accretion disk and channelled and collimated by the toroidal magnetic field 
built up by the rotation of the disk (Pudritz & Ray 2019). Since we now understand that a 
dynamically significant magnetic field can be present even in the earliest minihalos (Schober 
et al. 2012a, see also Section 3.3.2), an obvious question is whether the formation of Pop 
III stars is also associated with the presence of magnetized outflows. 

In practice, the answer to this question appears to depend on what one assumes re- 
garding the initial geometry of the field. Machida et al. (2008a), Machida & Doi (2013) 
and Sadanari et al. (2021) simulate the collapse of primordial gas clouds threaded by an 
initially uniform magnetic field and show that protostellar outflows are a common outcome 
unless the initial magnetic field strength is very small. On the other hand, the simulations 
performed by Sharda et al. (2020) and Prole et al. (2022a) that start with a tangled mag- 
netic field of the kind that we would expect to result from the turbulent dynamo find no 
evidence for the launching of jets or outflows. Whether this means that Pop III protostars 
do not produce outflows at all or simply that outflows are launched at a later point in the 
evolution of the system than is followed by these simulations remains to be determined. 

The other class of stellar outflows that are potentially of interest are winds driven by the 
Pop III stars themselves. At solar metallicity, massive stars can produce highly energetic 
winds that can have a profound impact on the surrounding gas (see e.g. Rahner et al. 2017). 
However, these winds are driven by the absorption of stellar photons by the lines of highly 
ionized metals and we expect the wind strength to be a function of metallicity (Kudritzki 
2002). They become weak or non-existent for zero metallicity (Krticka & Kubát 2006, 2009). 
Accounting for stellar rotation does not significantly change this picture: Stars rotating very 
close to their break-up velocity can drive strong winds, but lose so much angular momentum 
in the process that they do so only very briefly (Ekstróm et al. 2008). Stellar wind feedback 
therefore does not appear to be important for understanding the formation of Pop III stars. 


4.3. Large-scale impact 


In addition to the profound impact that feedback from the first stars has on the gas in 
their parent minihalos, feedback from these stars also strongly affects gas on much larger 
scales, both in the intergalactic medium and in neighboring protogalaxies. In this section, 
we review the most important of these impacts. 


www.annualreviews.org e The First Stars 


33 


34 


4.3.1. Photoionization. Ultraviolet photons that manage to escape from the halos hosting 
Pop III star formation will start to ionize the surrounding IGM. Prior to cosmic reionization, 
the Universe is highly opaque to these types of photons, and so this is essentially a local 
process: Individual star-forming protogalaxies will produce distinct ionized regions in the 
IGM, but no large-scale ionizing background will build up until after these individual regions 
begin to overlap, an event which does not occur until long after Pop II stars have taken 
over as the dominant source of ionizing photons (Hartwig et al. 2022). 

The impact of photoionization feedback from Pop III star formation on nearby minihalos 
and atomic cooling halos depends on their evolutionary state and on the duration of the 
feedback. Systems in which gas has already cooled and is in the process of forming stars 
are difficult to disrupt with external photoionization feedback, since the time taken for an 
ionization front to propagate through the system is much longer than the collapse time 
of the gas, once the gas density exceeds n ~ 100 cm? (see e.g. Whalen et al. 2008). On 
the other hand, systems that have not yet gathered together much cool gas can be strongly 
affected. Visbal et al. (2017) show for a constant ionizing flux that photoionization increases 
the minimum virial temperature required for collapse to Tu, ~ 40000 K if the ionizing flux 
exceeds Fart ~ 10° photons s -lcm- 2. For smaller values, photoionization is of limited 
importance, although cooling is suppressed in halos with Tu, < 10* K by the LW radiation 
that Visbal et al. (2017) assume accompanies the ionizing radiation (see Section 4.3.2 below). 
For Pop III minihalos with star formation efficiencies typical of those found in simulations, 
the ionizing flux exceeds Ferit only at distances of a few kpc or less, even when assuming 
an escape fraction of 100%. Moreover, the duration of Pop III star formation in individual 
minihalos is short, owing to the strong effects of photoionization and supernova feedback 
on small scales. Therefore, cooling will only be suppressed for a brief period of a few Myr. 
Over longer timescales, the effect of ionizing feedback within this region can actually be 
positive, since the residual fractional ioniztion in the recombining region will promote H> 
formation (O’Shea et al. 2005). Overall, we do not expect photoionization by Pop III- 
dominated sources to be the main source of feedback at high redshifts. Feedback from Pop 
II-dominated protogalaxies on surrounding minihalos is potentially more important, but a 
thorough discussion of this topic is outside of the scope of this review. 


4.3.2. Photodissociation. As discussed in Section 4.2.2, massive Pop III stars produce large 
numbers of LW band photons that can photodissociate H>. A large fraction of these escape 
from their parental minihalos (Kitayama et al. 2004; Schauer et al. 2015, 2017a), and they 
can propagate to large distances through the IGM owing to its very low opacity below 
13.6 eV (Haiman et al. 2000). The onset of Pop III star formation is therefore quickly 
followed by the build-up of an extragalactic LW background that pervades the Universe. 
This photodissociates Hz in newly-assembled minihalos, making it harder for the gas to 
cool (Haiman et al. 1997). Once the strength of the background exceeds a critical value, 
H» cooling is almost entirely suppressed, preventing gas in minihalos from collapsing and 
forming Pop III stars. 

The strength of the LW background required in order to suppress H> cooling in a 
minihalo can be estimated by a simple timescale argument (Oh & Haiman 2002). Models for 
the build-up of H> in primordial gas demonstrate that most molecules form within the first 
few recombination times, with the H> formation rate dropping substantially at later times 
owing to the loss of the required electrons and protons from the gas (Tegmark et al. 1997; Oh 
& Haiman 2002; Glover 2013). The LW background therefore only significantly suppresses 
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Hə formation once the photodissociation timescale, tais = 1/kais, becomes comparable to the 
recombination timescale, krec = 1/(krecn%o), where kais and krec are the photodissociation 
and radiative recombination rate coefficients, n is the characteristic number density of the 
gas and zo is its initial fractional ionization. For the conditions appropriate for a minihalo 
close to the minimum mass scale for collapse (see Section 3.1), that is T ~ 2000 K and 
zo ~ 2 x 1074, this yields a critical value (Oh & Haiman 2002) of 


(2) = 10 ^ fen, 24. 
n crit 


above which we expect the LW background to significantly affect the minihalo. Here, J21 is 
the strength of the background at the Lyman limit in units of 107?" ergs ^! Dg 


cm”? Hz”! sr}, 
n is the density of the gas at the center of the minihalo, and fz, is the self-shielding factor, 
ie. the ratio of the H» photodissociation rate to its value in the absence of self-shielding. 
Numerical simulations (e.g. Machacek et al. 2001; Yoshida et al. 2003; Kulkarni et al. 
2021) indicate that n ~ 10 cm”? is an appropriate value for the central density in the 
earliest minihalos, and so in the absence of self-shielding, we expect the LW background to 
become important once its strength exceeds a critical value of around J»1,«it ~ 1078. In 
this case, the LW background reduces the peak H> abundance by a factor proportional to 
Län deer OX Jal thereby increasing the minimum virial temperature required for efficient 
Ha cooling. At typical minihalo temperatures, the H> cooling rate scales approximately as 
An, X T^ and the cooling time scales as Lol cx (TŠ ras)”. Therefore, for the cooling time 
to remain approximately constant as Joi increases, we must have Tmin < J2/ 3, Moreover, 
since Merit ox Tu this implies that the minimum mass scale for collapse must increase as 
Merit < JE? for Joi > Jan sain, 

In the upper panel in Figure 10, we show a compilation of predictions for how Merit varies 
as a function of Jz; taken from a number of studies in the literature. These studies range 
from simple semi-analytical models (‘Trenti & Stiavelli 2009; Visbal et al. 2014a) to detailed 
numerical simulations both with (Skinner & Wise 2020; Kulkarni et al. 2021; Schauer et al. 
2021) and without (Machacek et al. 2001; Wise & Abel 2007b; O'Shea & Norman 2008) 
H» self-shielding. Although there is considerable scatter,? we see that the semi-analytical 
models and the simulations without H> self-shielding yield results that broadly agree with 
our expectations based on the discussion above. However, the situation becomes more 
complicated once we account for Hə self-shielding. Schauer et al. (2021) include this in 
their simulations, but do not find it to have a large effect. They report values of Merit 
that are a factor of 2 — 3 smaller than those predicted by simulations that do not account 
for self-shielding, but whether this difference is due to the inclusion of self-shielding or due 
to other differences in the chemical or thermal treatment is unclear. On the other hand, 
Skinner & Wise (2020) and Kulkarni et al. (2021) find He self-shielding to have a profound 
impact on Merit, keeping it below 107 Mo for LW background strengths of up to J21 ~ 10. 
If true, these results would imply that the LW background does not play an important role 
in regulating high redshift star formation, since the mean strength of the background does 
not exceed J21 = 1 until after Pop II star formation has already started to dominate (see 
e.g. Trenti & Stiavelli 2009; Ahn et al. 2009). 


?'This scatter likely results from a combination of differences in the microphysical treatment (e.g. 
selection of reactions, choice of reaction rate coefficients and cooling function) and differences in 
the definition of Mcrit in the different studies. 
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Figure 10 


Top: Minimum minihalo mass required for Pop III star formation (Merit) plotted as a function of 
J21, the strength of the LW background in units of 107?! ergs-! cm~? Hz”! sr—!. The values 
shown are taken from studies by Machacek et al. (2001) (MBAOI; blue plus signs), Wise & Abel 
(2007b) (WA07; black crosses), O’Shea & Norman (2008) (ON08; red circles), Trenti & Stiavelli 
(2009) (TS09; blue line), Visbal et al. (2014a) (black points), Skinner & Wise (2020) (SW20; blue 
pentagons), Kulkarni et al. (2021) (KVB21; green squares) and Schauer et al. (2021) (green 
triangle). For TS09, we show the values computed assuming z = 20. SW20 find H2 cooling and 
star formation occurring in minihalos with a broad range of masses for 107! < Jo1 < 1; we show 
only the few least massive points from their study here. Finally, for S21, we show the values they 
report for the lowest mass at which minihalos form stars. The average mass at which minihalos 
first form stars (their Mave) is about a factor of three larger than this minimum value. Bottom: 
Evolution of Merit with redshift for a selection of models that account for the combined effects of 
baryonic streaming and the LW background. The lines are computed for a streaming velocity 

Ust = 0.8 (in units of the rms value) and the shading indicates the impact of varying this by +0.5. 
For the time evolution of the LW background, we adopt a simple fit to the values computed by 
Greif & Bromm (2006). The meaning of the various lines is the same as in Figure 1. 


Another view of the current uncertainty in the impact of LW feedback is given in the 
bottom panel of Figure 10, where we show the redshift evolution of Merit in several different 
models, assuming a streaming velocity vs, = 0.8 (in units of the root-mean-squared value) 
and a LW background that evolves as Jo1 = 10?-7/ ZA which is a reasonable fit at z > 10 
to the values computed by Greif & Bromm (2006). At z > 25, the LW background is 
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very weak and the evolution of Merit in all of the models is essentially the same as in the 
Ju = 0 case discussed earlier (see Section 3.1.2 and Figure 1). At z < 25, however, the 
growing background drives a substantial increase in Merit in the Schauer et al. (2021) and 
F13 models, but has little impact on it in the Kulkarni et al. (2021) and OHS models. 

Figure 10 also demonstrates that in halos more massive than ~ 107 Mo at z = 20, Merit 
becomes highly insensitive to Jan even in the absence of effective H> self-shielding. The 
reason for this change of behavior is that these halos have high enough virial temperatures 
to allow them to cool initially via Lyman-o cooling. In these atomic cooling halos, the 
gas can reach much higher densities without the need for H> cooling, meaning that a much 
stronger LW background is needed to completely suppress H> formation. The key density in 
this case is Nerit ~ 104 cm °, the critical density above which the rotational and vibrational 
level populations of H> start to approach their local thermodynamical equilibrium (LTE) 
values. Above nerit, the H> cooling time becomes independent of density, while the free- 
fall collapse time continues to decrease with increasing density. In addition, collisional 
dissociation of H> also becomes more effective in gas with n > nci. Together, these facts 
lead to a bifurcation in the behavior of atomic cooling halos depending on the amount of 
H> they have managed to form by the time their density reaches nerit. If they have formed 
enough H> to provide effective cooling, then the gas remains cold as it collapses, with 
its subsequent evolution being very much the same as in the standard Pop III formation 
scenario (Section 3.2). On the other hand, if H> cooling has not become effective by the 
time the gas reaches Merit, it will never become so at higher densities. Instead, the gas 
will remain warm, with cooling arising primarily from bound-bound transitions in atomic 
hydrogen and from H^ formation (Omukai 2001; Schleicher et al. 2010). 

'The combined effect of the higher temperatures and higher characteristic densities en- 
countered in atomic cooling halos is a dramatic increase in the value of Jan oc. One-zone 
models of the thermal and chemical evolution of highly irradiated gas in atomic cooling 
halos typically find J21,crit ~ 1000 if the spectrum of the LW background is dominated by 
Pop III stars (see e.g. Sugimura et al. 2014; Glover 2015a,b). Three-dimensional models 
find even larger values, J21,crit ~ 10* (Shang et al. 2010; Latif et al. 2015; Hartwig et al. 
2015), although the scatter between different simulations can be considerable. 

'To put these numbers into context, we expect the mean strength of the LW background 
at reionization to be Joi ~ 1 (Haiman et al. 2000; Ahn et al. 2009). Therefore, the LW 
background is unable to suppress H> cooling in the vast majority of atomic cooling halos, 
although it does produce systematically warmer temperatures and higher accretion rates 
than would be the case without the background (Latif & Volonteri 2015; Regan & Downes 
2018). Complete suppression of H> cooling occurs only in a small number of atomic cooling 
halos that happen to lie very close to a strong source of LW photons (Dijkstra et al. 2008; 
Agarwal et al. 2012; Dijkstra et al. 2014; Visbal et al. 2014b; Regan et al. 2017; Agarwal 
et al. 2019) or that are undergoing strong dynamical heating from a succession of mergers 
(Regan et al. 2020). However, despite their rarity, these systems are of considerable interest 
as the possible formation sites of supermassive Pop III stars (see Section 3.4). 


4.3.3. X-rays. Although Pop III stars themselves do not produce significant emission at 
X-ray wavelengths, we nevertheless expect Pop III star formation to be closely associated 
with the production of X-rays. We discuss the origin of this emission as well as its impact 
onto the surrounding gas. 
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Stellar sources of high-redshift X-rays. Three main mechanisms potentially contribute 
to the X-ray emission from Pop III stars. First, if some fraction of these stars are located in 
tight binaries, then these systems can become high-mass X-ray binaries (HMXBs, see e.g. 
Mirabel et al. 2011). Second, hot shocked gas in Pop HI supernova remnants (SNRs) emits 
strongly at X-ray energies until it cools. Finally, Pop III SNRs may also produce X-rays 
indirectly, by acting as a source of cosmic ray electrons. Because of the high CMB energy 
density at high redshift, cosmic ray electrons with energies of up to a few hundred GeV lose 
most of their energy by inverse Compton scattering CMB photons, producing a power-law 
inverse Compton spectrum that extends up to X-ray energies (Oh 2001). 

A useful metric for comparing the relative importance of these mechanisms is the X-ray 
luminosity they produce per unit star formation rate, Lx /SFR. In the local Universe, at 
metallicities close to solar, HMXBs produce a luminosity in the 0.5 — 2keV energy band 
given by Lx /SFR z 6 x 10? ergs”! (Maar 1) 1 (see e.g. Glover & Brand 2003; Kaur et al. 
2022), although there is roughly a factor of two scatter between different observational 
determinations. However, the contribution from SNRs‘ is less well determined. Estimates 
of the ratio of the X-ray luminosity to the mechanical luminosity range from values as small 
as 2 x 1074 (e.g. Helfand & Moran 2001) to ones as large as a few percent (e.g. Smith 
et al. 2005; Franeck et al. 2022), albeit with substantial time variability and scatter from 
region to region. If we take a value near the midpoint of this range, Lx/Lmech = 107? as 
a representative average (Rogers & Pittard 2014; Franeck et al. 2022), then we find that 
Lx/SFR z 6 x 1095 ergs ! (Mayr!) 1 for SNRs (Leitherer et al. 1999), i.e. around 10% of 
the HMXB contribution. Therefore, on average in the local Universe, X-rays from HMXBs 
dominate over those produced by hot gas in SNRs. Accounting for the inverse Compton 
contribution does not substantially change this picture (Oh 2001). 

How do these conclusions change as we move to higher redshift? Both models and 
observations show that HMXBs produce X-rays more efficiently as the metallicity drops, 
with Lx/SFR increasing by an order of magnitude or more as the metallicity decreases 
for 1 Zo to 0.03Ze (Kaur et al. 2022). At even lower metallicities, we have no direct 
observational constraints, and theoretical models disagree on whether Lx /SFR asymptotes 
to a constant value (Fragos et al. 2013) or continues to increase (Brorby et al. 2016). It 
is also unclear whether this finding of a higher Lx/SFR remains valid if we extrapolate 
all the way down to Z — 0, with results from one study indicating that HMXBs formed 
from Pop III stars may actually be less effective at producing X-rays than in the solar 
metallicity case (Sartorio et al. 2023), although we note that this result depends sensitively 
on the choice of the Pop III IMF. Turning our attention to SNRs, we note that we do not 
expect the fraction of the energy they radiate as X-rays to be particularly sensitive to the 
metallicity. However, if the Pop III IMF is dominated by massive stars (see Section 3.2.3), 
then the number of SNe that explode per unit star formation rate will be a factor of a 
few larger than in the local Universe. In addition, if a significant number of Pop III stars 
explode as hypernovae or pair-instability supernovae (see Section 4.3.4 below), then the 
energy available per explosion will also be larger, implying a higher X-ray luminosity even 


^We do not consider stellar winds here, as even at solar metallicity, these do not succeed in 
producing large X-ray luminosities (Wang & Helfand 1991; Harper-Clark & Murray 2009), due to 
some combination of the escape of hot gas from wind-driven bubbles and efficient cooling of the 
gas at the bubble boundary (Lancaster et al. 2021). At low or zero metallicity, stellar winds will be 
much weaker and hence will be even less effective at producing X-rays. 
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if Lx / Lech remains constant. In addition, the inverse Compton contribution is potentially 
far more important here than at low redshift. For a typical supernova, around 0.1 — 0.2% 
of its initial kinetic energy is used to accelerate relativistic electrons (Lacki & Thompson 
2013) and so if the bulk of this is converted to X-rays by inverse Compton scattering, the 
result is a contribution to Lx /SFR similar in size to that from hot gas. 

In summary, if Pop III HMXBs are as efficient at producing X-rays as present-day 
HMXBs and if most Pop III SNe are standard core-collapse supernovae with explosion 
energies Esn ~ 10?! erg, then the majority of the X-ray emission associated with Pop III 
star formation will come from HMXBs. For this reason, it is this scenario that has attracted 
most attention in the literature, as we discuss in more detail below. However, scenarios in 
which Pop III supernovae dominate the production of X-rays at high redshift cannot be 
completely ruled out (see e.g. Ricotti 2016, for an example). 


Impact of the high-redshift X-ray background. The X-rays produced by Pop III 
HMXBs and other sources can propagate to large distances through the IGM. The Universe 
is optically thin to X-ray photons with energies E > 1.8[(1 + z)/15] /? 
ing photon mean free path exceeds 1 Mpc for photon energies E > 0.25 keV (Furlanetto 
et al. 2006). Therefore, once Pop III star formation begins, a pervasive X-ray background 


keV and the comov- 


develops. It affects gas in high redshift minihalos in two ways. First, it partially ionizes 
it, thereby catalyzing the formation of additional H> (Haiman et al. 2000). Second, it 
also heats the gas. The first of these effects acts as positive feedback, making cooling and 
star formation more likely, while the second one acts as negative feedback. Their relative 
importance depends on the density of the gas and the strength and spectral shape of the 
background. Heating dominates at low densities and can completely suppress the collapse 
of gas into low-mass minihalos if it is able to heat the IGM to a temperature above the virial 
temperature of the minihalo. The strength of the background required to accomplish this 
depends sensitively on the adopted low-energy cutoff, since the heating rate due to photons 
of energy E scales approximately as E^? owing to the strong energy dependence of the 
photionization cross-section. For example, Hummel et al. (2015) take a minimum photon 


energy of 1 keV and estimate that a critical value of Jx,crit © 10 21 ergs 1cm”? Haler) 


at 1 keV is required to suppress the collapse of gas into a 10° Mç halo. On the other hand, 
Park et al. (2021a) adopt a minimum photon energy of 200 eV and find that a much weaker 
background is sufficient to suppress collapse, Jx,crit 10 28 ergs 1cm”? Hz”! sr”!. 


At higher densities, the heating provided by the K-rays becomes unimportant in com- 


parison to the increased cooling provided by the additional Hz. Consequently, the gas cools 
to lower temperatures than in the case without X-rays (Machacek et al. 2003; Hummel 
et al. 2015; Park et al. 2021a). The density at which the behavior changes lies in the range 
n ~ 1— 100 cm `Š, with the value depending on the strength of the background. However, if 
the X-ray background is very strong, there is so much heating at low densities that the gas 
never reaches the positive feedback regime, and the feedback in this case is wholly negative. 

At even higher densities, the gas becomes optically thick to the X-ray background, which 
therefore plays no direct role on scales comparable to the size of the protostellar accretion 
disk. However, its presence indirectly affects the behavior of the gas on these scales, as the 
lower temperatures that are reached on larger scales lead to a slower inflow of mass to the 
center of the halo. This may result in the formation of a less massive, more stable accretion 
disk (Hummel et al. 2015; Park et al. 2021a,b), which means fewer fragments and implies a 
more top-heavy IMF as well as lower stellar multiplicity. 
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Another interesting question is whether or not positive feedback from X-rays dominates 
over negative feedback from LW photons. The answer obviously depends on the relative 
strength of the two backgrounds (Park et al. 2021a), but models that self-consistently follow 
their build-up typically find that LW feedback dominates (see e.g. Glover & Brand 2003). 

Finally, the heating of the IGM produced by even a weak X-ray background has an 
important impact on its visibility in the 21 cm line of atomic hydrogen (Furlanetto et al. 
2006). Prior to the onset of Pop III star formation, the IGM is colder than the CMB and 
hence the 21 cm signal is visible in absorption (see Section 5.1.3). However, once X-rays heat 
the IGM above the CMB temperature, the absorption signal vanishes and the line instead 
becomes visible in emission. Measuring the redshift at which the transition from absorption 
to emission occurs therefore allows us to constrain the strength of the high redshift X-ray 
background, and may also allow us to distinguish between the contributions made to it by 
Pop III and Pop II sources (Mirocha et al. 2018). 


4.3.4. Supernovae. In this section, we discuss how the final fate of a Pop III star depends 
on its initial mass, with a particular focus on identifying the range of masses that lead to 
supernova explosions. We also briefly outline the impact that these explosions have on the 
surrounding gas. 


The fate of massive Pop III stars. How massive Pop III stars end their lives depends 
on the mass of their helium core at the moment that the star dies (Heger et al. 2003). For 
non-rotating stars with negligible mass loss, there is a simple mapping between helium core 
mass and initial mass, and so it is common to discuss the fates of these stars as a function 
of their initial mass. Below M ~ 9Mo, Pop III stars end their lives as white dwarfs, 
supported by electron degeneracy pressure. Above this mass, electron degeneracy pressure 
is not sufficient to support the star. Stars in the narrow mass range 9Mo < M < 10 Mo 
form degenerate oxygen-neon cores that collapse, while those with 10 Mç < M < 100 Mo 
form iron cores that also undergo collapse. In stars with 10Mo < M < 25 Mo, this core 
collapse triggers a supernova explosion (hereafter referred to as a core-collapse supernova 
or CCSN) with an energy of around 10?! erg that leaves behind a neutron star remnant. 
For 25 Mo < M < 40 Mo, a core-collapse explosion also occurs, but is unable to completely 
unbind the stellar envelope from the neutron star remnant. Instead, the inner portions 
of the envelope fall back onto the neutron star, causing it to collapse and form a black 
hole. The energy released into the surrounding environment depends on the extent of this 
“fallback” and can potentially be much smaller than 10°! erg, in which case we refer to 
the CCSN in question as a faint supernova. For masses between 40 Mọ and a somewhat 
uncertain upper limit in the range 70 — 100Mo, no supernova explosion occurs. Instead, 
the star simply collapses directly into a black hole at the end of its life. 

For more massive stars, a new physical effect becomes important. These stars have 
such high central temperatures during the later stages of their evolution that they become 
susceptible to the so-called pair instability (Bond et al. 1984; Heger & Woosley 2002). 
Briefly, as the temperature increases, so too does the fraction of photons in the core that 
have energies greater than 1.022 MeV, the rest mass of an electron-positron pair. Interaction 
of these photons with the surrounding nuclei can lead to their conversion into electron- 
positron pairs. Since the core is supported by radiation pressure, the loss of photons via 
this mechanism results in a reduction in the pressure and hence contraction of the core. 
This further increases the core temperature, making pair production more likely, and in 
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sufficiently massive stars, this can become a runaway process triggering the dynamical 
implosion of the core (Woosley 2017). This implosion is terminated by the onset of rapid O 
and Si burning. Early calculations predicted an onset of the pair instability above a mass 
of around 100 Mo (Heger et al. 2003), but more recent calculations find lower values of 
~ 70 Mç (Woosley 2017), albeit with some remaining uncertainty due to the uncertainty 
in the rate of the !?C(o, y)!9O nuclear reaction (Farmer et al. 2019). In stars above this 
mass, but below 140 Mo, this leads to a pulsational instability in which the star undergoes 
a series of implosions, nuclear flashes and re-expansions that disrupt the stellar envelope 
and part of the core. This process terminates once the core mass becomes too small for it 
to be susceptible to the pair instability, following which it evolves in the same fashion as 
a lower mass star. The energy released by this instability ranges from < 10?! erg to a few 
times 10?! erg, and the ejected mass also spans a large range of values (Woosley 2017). 

Above 140 Mo, pulsations do not occur, as the energy released in the initial nuclear 
flash completely disrupts the star, leaving no remnant behind. The resulting explosion is 
termed a pair-instability supernova (PISN). It has a typical energy of around 10”? erg, two 
orders of magnitude larger than the value of a typical CCSN. Because all of the mass of the 
star is expelled in the ejecta, PISNe typically have very large metal yields. Finally, for Pop 
III stars more massive than 260 Mo, the pair instability occurs, but the energy released 
by the resulting rapid nuclear burning is not sufficient to disrupt the star, which instead 
simply collapses to a black hole (Heger & Woosley 2002). 

Accounting for the effects of rotation has two main effects on the picture sketched 
above. First, rotation-driven mixing leads to rapidly rotating Pop III stars having larger 
helium core masses than non-rotating stars of the same mass (Yoon et al. 2012). Since 
the onset of the pair instability depends on the core mass, an important consequence of 
this is a reduction in the minimum mass required for a PISN, which can reach values as 
low as M ~ 100 Mo for rapidly rotating stars. Second, rapidly rotating Pop III stars can 
potentially explode as jet-driven supernovae (see e.g. Grimmett et al. 2021, and references 
therein). These are widely believed to be responsible for the unusually bright class of SNe 
known as *hypernovae" (HNe, see Umeda & Nomoto 2002), which have typical energies of 
10?? erg. Rapid rotation is also a necessary (but not sufficient) condition for a supernova 
to produce a gamma-ray burst (GRB, see Toma et al. 2016). 

Unfortunately, the rotational properties of Pop III stars are poorly constrained. In 
order to measure the angular momentum of Pop III protostars reliably in simulations, 
extremely high resolution is required in order to distinguish between tight binaries with 
separations of a few AU and single stars. One of the few sets of simulations that do have 
sufficient resolution are those presented by Greif et al. (2012). Stacy et al. (2013) analyzed 
the rotational properties of the protostars formed in these simulations and showed that 
most were rotating very rapidly, with rotational velocities between ~ 5096 and 100% of the 
Keplerian velocity. However, these simulations cover only a very short period (~ 10 yr) at 
the start of the evolution of the protostellar system and do not account for the effects of the 
magnetic field, which may provide strong magnetic breaking on protostellar scales (Hirano 
& Bromm 2018). In view of these uncertainties, most studies of the impact of Pop III SNe 
on their surroundings have focussed on the non-rotating case. 


Impact of Pop III supernovae. Pop III supernovae affect their surroundings via both 
mechanical feedback (the injection of energy and momentum) and chemical feedback (the 
enrichment of the gas with metals). We first consider the role played by mechnical feedback. 
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The gravitational binding energy of a minihalo with a mass close to Merit is ~ 10°! erg, 
comparable to the energy of a single CCSN. We would therefore expect a single CCSN to 
be able to eject a considerable fraction of the gas from a low-mass minihalo, provided that 
the supernova remnant does not lose too much energy via radiative cooling. 

In practice, the main factor that determines whether or not radiative cooling of the 
supernova remnant is effective is the density of the region in which the supernova explodes. 
If prior photoionization feedback has cleared away most of the gas, then cooling is ineffective 
and even a single supernova can eject a large fraction of the gas from the halo (Kitayama & 
Yoshida 2005; Ritter et al. 2012; Jeon et al. 2014; Chiaki & Wise 2019). On the other hand, if 
prior photoionization feedback is ineffective, as would be the case for e.g. a 10—15 Mç CCSN, 
then little gas is ejected (Jeon et al. 2014; Chiaki & Wise 2019). Highly energetic SNe, such 
as hypernovae or PISNe, are less sensitive to the local environment and typically clear 
almost all gas out of the halo (Greif et al. 2010; Magg et al. 2022c). This picture changes in 
massive minihalos, due to their larger gravitational binding energies. For minihalo masses 
larger than a few times 10? Mo, a single CCSN is no longer sufficient to eject the gas, 
although a single PISN remains effective (Chiaki & Wise 2019; Magg et al. 2022c). 

The time taken for gas to become available again for star formation after the explosion 
of a Pop III SN, the so-called fallback or recovery time, varies considerably with halo mass 
and SN progenitor mass, ranging from a few Myr to 100 Myr or more. Smaller values 
are typically associated with lower mass progenitors or higher mass minihalos, but there is 
considerable scatter from halo to halo even for fixed progenitor mass. Magg et al. (2022c) 
argues that this is a consequence of the small scale structure of the gas prior to the supernova 
explosion. The greater the amount of gas located in dense clumps, the better it is able to 
survive the impact of the SN remnant and the shorter the resulting fallback time. 

Supernova remnants that successfully expel gas from their parent minihalo will sub- 
sequently expand into the surrounding IGM, reaching final sizes of 500 — 1500 pc by the 
time that they finish expanding (Greif et al. 2010; Jaacks et al. 2018), depending on their 
initial energy and the density of the surrounding gas. These sizes are large enough to reach 
neighboring minihalos. If these minihalos are already in the process of forming stars, then 
the interaction with the expanding supernova remnant will have little effect on them. It will 
strip away the low density gas in the outskirts of the minihalo, but will not penetrate into 
the dense core (Chen et al. 2017). Minihalos that are still assembling their gas are more 
strongly affected, as in this case the expanding remnant can strip most of the gas from 
the halo, preventing it from cooling and collapsing (Sigward et al. 2005). However, these 
minihalos will already have been irradiated by photoionizing radiation from the supernova 
progenitor (Jaacks et al. 2018), which has a similar suppressive effect, and so this direct 
mechanical feedback is generally considered to be of limited importance. 

Turning to chemical feedback, we note that this proceeds in two qualitatively different 
ways: internal enrichment (i.e. enrichment of gas in the halo where the supernova exploded) 
and external enrichment (i.e. enrichment of nearby halos that interact with the supernova 
remnant). As we might expect from the discussion above, the relative importance of these 
mechanisms depends on how effectively supernovae can drive gas and metals out of their host 
halos, with internal enrichment dominating in more massive halos and external enrichment 
dominating in lower-mass halos (Chiaki & Wise 2019). The likelihood of external enrichment 
occurring also depends on the local number density of minihalos capable of forming stars, 
as this determines how many are likely to be located within the region of the IGM enriched 
by the supernova. Therefore, we would expect the importance of external enrichment to 
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also depend on Merit, with this process becoming less important as Merit increases. As an 
illustration, consider the studies of early metal enrichment by Hicks et al. (2021) and Magg 
et al. (2022c). Hicks et al. (2021) find that ~ 60% of the metal-enriched minihalos in their 
simulation gain their metals via external enrichment. However, these halos are small, with 
most having masses M < 10°Mo, and it is unclear how many will actually form stars. 
In comparison, Magg et al. (2022c) examine only the enrichment of halos that successfully 
form stars and find that in this case, only around 3596 are externally enriched. Although 
these studies are a good start, more work along these lines is necessary to fully pin down 
the relative importance of external versus internal enrichment and how this depends on the 
local extragalactic environment, the strength of the LW background etc. 

Finally, a quantity of great interest for the interpretation of observations of extremely 
metal-poor stars in our own galaxy (see Section 5.2.1) is the metallicity of the gas enriched 
by a single Pop III supernova. This depends on the efficiency with which metals are mixed 
into the surrounding environment: the greater the amount of mixing, the lower the resulting 
metallicity. However, we can put an upper limit on this metallicity by assuming that metals 
only mix within the volume occupied by the supernova remnant at the time that it finishes 
expanding. Magg et al. (2020) show that in this case, the minimum mass of gas into which 
the metals are mixed is given by 


Maitmin = 1.9 x 10* ER n Mo , 25. 


where Eau is the kinetic energy of the SN explosion in units of 10?! erg and n is the number 
density of the gas in which the explosion takes place in units of cm~?. For explosions taking 


—011 4 | and so a CCSN will mix its metals into 


place within a pre-existing HII region, n 
a minimum of 2 x 10! Mç of gas, while a PISN will mix its metals into a minimum of 
2x 10° Mo of gas. The total mass of metals ejected by a Pop III CCSN ranges from around 
0.1 to a few Mc (Nomoto et al. 2006), and taking a value at the top end of this range 
then yields a maximum metallicity of Zmax ~ 0.01 Zo for enrichment by a single CCSN. 
For a PISN, the mass of metals ejected is much larger, ~ 100 Mc (Heger & Woosley 2002), 
but since this is mixed into a much larger gas mass, the resulting maximum metallicity 
is actually smaller, Zmax ~ 3 x 10? Zo. Magg et al. (2020, 2022c) show that in practice, 
individual Pop III supernovae mix their metals into gas masses ranging from a few times this 
minimum value to as much as 10* times more. Large values typically correspond to cases of 
external enrichment, when a large mass of pristine gas associated with an existing minihalo 
is enriched by a small mass of metals from an external source. It should also be stressed 
that this mixing is generally not uniform, there can be substantial metallicity variations 
within the enriched volume (Ritter et al. 2015; Tarumi et al. 2020). Individual Pop III 
SNe can therefore produce metal-enriched gas with a broad range of metallicities, making 
it plausible that at least some of the extremely metal-poor stars observed in the Milky Way 
will have formed from gas enriched by only a single Pop III SN (see also Section 5.2.1). 


4.4. Transition to Pop II star formation 


Second generation stars, sometimes termed early Pop II stars, formed from material that 
was enriched by metals from the first generation of stars. Unlike Pop III stars, which 
have not yet been directly detected (Section 5.1.1), members of the second generation have 
already been found in surveys looking for extremely metal-poor stars in our Milky Way and 
neighboring satellite galaxies (Section 5.2.1). As discussed in Section 3.2.1, the chemical 
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composition of the gas governs its thermodynamic response and therefore its fragmentation 
behavior under different environmental conditions. In turn, this determines the stellar mass 
spectrum (Section 3.2.2), which for genuine Pop III stars is predicted to be very different 
compared to what we find in the Universe today (Section 3.2.3). Unlike purely primordial 
gas, where cooling is largely provided by hydrogen (see top right of Figure 3), metal- 
enriched gas has access to a much wider range of cooling mechanisms, allowing it to reach 
lower temperatures at various points in its evolution. This enhances the fragmentation of 
the gas and reduces the rate at which the fragments accrete additional material. Together, 
these effects are widely believed to lead to an IMF dominated by low-mass stars, similar to 
the present-day IMF (Kroupa 2002; Chabrier 2003). However, there remains the question 
of when and where this transition occurs, whether it is sudden or gradual, and what is 
the underlying physical cause for it. Indeed, there are two competing models based on two 
different low-metallicity cooling channels, each associated with a different critical metallicity. 

The first model is based on atomic fine-structure lines from alpha elements such as 
carbon or oxygen. These start to dominate the coolingat number densities around 10* cm”? 
and above a critical metallicity of about 107? Ze (see, e.g. Chon et al. 2021). This leads 
to a markedly different evolution in the log T—log n diagram as indicated in the upper 
right part of Figure 3. The required metallicity (or, more precisely, the required carbon and 
oxygen abundances) can be estimated by comparing the cooling rate due to metals with the 
compressional heating rate at a characteristic temperature of 200 K and density of 10* cm”? 
(see e.g. Bromm & Loeb 2003; Santoro & Shull 2006), corresponding to the point at which H> 
cooling becomes inefficient in the Pop III case (see Figure 3). To quantify this, Frebel et al. 
(2007) introduced the transition discriminant, Dtrans = logy DER +0.3 x wn and 
argued that low-mass stars can only form for metallicities Dtrans 2 —3.5. 

The second line of reasoning considers dust cooling as the primary agent for fragmenta- 
tion, which leads to a considerably lower critical metallicity in the range 10^? to 1079 Zo. 
This value was first estimated using analytical or simple one-zone numerical models that 
compared the strength of various cooling and heating processes in a metal-poor environment 
(e.g. Omukai et al. 2005; Schneider et al. 2006, 2012a; Chiaki et al. 2013b,a). However, it 
has subsequently received strong support from detailed 3D numerical simulations of dust 
cooling in dense, metal-poor gas (e.g. Tsuribe & Omukai 2006, 2008; Clark et al. 2008; 
Dopcke et al. 2013; Safranek-Shrader et al. 2014; Chiaki et al. 2016; Chiaki & Wise 2019; 
Chon et al. 2021; Chiaki & Yoshida 2022). Together, these models suggest that dust is the 
driving agent for the transition in the IMF. It occurs at a metallicity of about 10^? solar 
and leads to an IMF peaking below 1 Mo with a functional form similar to what is found 
in the solar neighborhood (Kroupa 2002; Chabrier 2003). 

The existence of the star SDSSJ1029151+4172927 (Caffau et al. 2011a) provides strong 
evidence for the validity of the dust-induced fragmentation model (e.g. Schneider et al. 
2012b; Klessen et al. 2012; Chiaki et al. 2014; Bovino et al. 2016) and is a challenge to 
models based on fine-structure cooling. SDSS J10291514-172927 is a truly primitive star in 
the constellation of Leo with a metallicity that falls well below the suggested critical value 
of Dtrans for atomic fine-structure lines to be relevant. The star has elemental abundances 
in the range 10~° to 1074 of the solar value for all of the elements measured in its spectrum 
(Caffau et al. 2012), setting it apart from other extremely metal-poor stars which typically 
have enhanced CNO abundances despite being very iron poor (Beers & Christlieb 2005). 

What is the physical origin of this behavior? Simple numerical experiments (Li et al. 
2003) indicate that turbulent gas in a cooling regime (i.e. one in which T' drops with in- 
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Figure 11 


Illustration of the change of fragmentation behavior and resulting stellar mass spectrum for 
different metallicities, ranging from 10-9 Z to 107! Ze. For high metallicities the resulting IMF 
is bottom heavy with a peak below one Mọ and a power-law decline above, similar to the 
present-day IMF. At low Z the mass spectrum becomes logarithmically flat and top heavy similar 
to what we expect for the Pop III IMF (cf. Figure 6). This is a reproduction of Figure 10 of Chon 
et al. (2021) and used here with permission of the authors and MNRAS. 


creasing density and the effective adiabatic index is eg < 1) can fragment very efficiently 
to build up a cluster of low mass stars. Conversely, star formation in a heating regime (T' 
increasing with density, Yer > 1) is biased towards forming very few and very massive ob- 
jects. It has been posited (Larson 2005) that the transition between these regimes provides 
ideal conditions for efficient fragmentation and introduces a characteristic mass scale to the 
system. This mass scale is approximately the Jeans mass at the density and temperature 
of the transition point (Jappsen et al. 2005). It can be understood as a consequence of the 
fact that fragmentation in turbulent, self-gravitating gas often occurs in filaments, which 
are unstable against fragmentation for yeg < 1 and stable for ye > 1 (Larson 1985), with 
isothermal gas constituting a critical case (Kawachi & Hanawa 1998). As all current models 
of the assembly of the first galaxies indicate that the gas develops a filamentary morphology, 
it is reasonable to assume that filament fragmentation will dominate in the earliest galaxies 
just as it does in the local ISM of the Milky Way (Hacar et al. 2022). 

The transition from yer < 1 to yee > 1 corresponds to a downward “kink” in the effective 
equation of state. Inspection of Figure 3 reveals two sets of such kinks, one associated 
with atomic fine-structure line cooling (for metallicities —3 < logio Z/Zc < 0) and one 
with dust cooling (in the range —5 < logio Z/Zo < 0). The associated characteristic 
masses are ~ 100 Mo for fine structure cooling and « 1 Mo for dust cooling. This implies 
that atomic fine-structure cooling is not able to produce the extremely metal-poor low- 
mass stars we detect in stellar archeological surveys (see Section 5.2) in sufficiently large 
numbers, providing additional evidence for dust cooling being the dominant process driving 
the transition from the Pop III to the Pop II IMF 

'To summarize the above, in Figure 11 we show the stellar mass spectrum that results 
from systematically varying the metallicity of the star-forming gas (for further details, see 
Chon et al. 2021). Although similar caveats apply here as for our previous discussion of the 
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Pop III IMF, we nevertheless see a clear influence of metallicity. Values close to solar lead to 
a bottom-heavy IMF which peaks around 0.2—0.3 Mo with a power-law drop towards larger 
masses, similar to the present-day values (Kroupa 2002; Chabrier 2003). As the metallicity 
gets smaller, deviations from the present-day IMF become more pronounced and the mass 
spectrum becomes logarithmically flat and top heavy, more comparable to the Pop HI IMF. 


5. Observational probes of Pop III stars 


Characterizing the properties of Pop III stars is difficult, as most of what we know about 
them comes from theoretical models and numerical simulations rather than observational 
constraints, which are hard to obtain and highly indirect at best. In this section, we provide 
an overview of the available data from the high-redshift and the present-day Universe. 


5.1. High-redshift Universe 


There are several different approaches that we can take when trying to use observations of 
the high-redshift Universe to constrain the properties of Pop III stars. First, we can try to 
detect the stars themselves. Second, we can observe the bright explosive transients that the 
most massive Pop III stars produce at the end of their lives (supernovae and gamma-ray 
bursts), Finally, we can look for larger-scale signatures of the energetic feedback produced 
by massive Pop III stars, particularly through the imprint they leave on the cosmic 21-cm 
background. 


5.1.1. Direct detection of Pop III stars. Many of our uncertainties regarding the formation 
of Pop III stars could be resolved if we could observe them directly in the high-redshift 
Universe. The possibility of doing so with current or near-future facilities has been investi- 
gated by a number of authors (see e.g. Oh et al. 2001; Scannapieco et al. 2003; Greif et al. 
2009; Zackrisson et al. 2011, 2012; Rydberg et al. 2013; Mas-Ribas et al. 2016; Riaz et al. 
2022). Unfortunately, the general conclusion of these studies is that this is not possible, at 
least if we are interested in Pop III stars with masses < 1000 Mo forming in minihalos or 
atomic cooling halos with efficiencies similar to those predicted by numerical simulations. 
A recent calculation by Schauer et al. (2020) highlights the problem. They show that a 
single 1000 Mo Pop III star forming at a redshift z > 10 will produce a flux in the NIR- 
CAM wide-field filters on JWST of between 5 x 107“ and 2 x 107? nJy, considering only 
the stellar emission. Accounting for nebular emission increases these numbers by about a 
factor of ten, but they are still far below the limiting values of a few nJy reachable with 
deep JWST imaging (see e.g. Merlin et al. 2022), although they would be in reach of a 
hypothetical 100 m extraterrestrial telescope (Schauer et al. 2020). 

In order to be detectable by JWST, our source therefore needs to be about 100 times 
brighter than we expect for a single very massive Pop III star. One way to reach the 
necessary brightness would be to have a cluster of such stars. However, the required stellar 
mass — around 10? Mọ in massive Pop III stars (Zackrisson et al. 2011) — does not seem 
achievable in a minihalo, given what we know about Pop III star formation. For instance, 
simply building up this much mass in a central cluster given the typical mass inflow rates we 
measure in simulations would require ~ 10 Myr, but feedback from the assembling cluster 
will disrupt the inflow long before this. Moreover, the required star formation efficiencies 
are unreasonably large, for instance, a minihalo with a mass close to Merit would have to 
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convert close to 100% of its baryonic mass into stars. Clusters of the required size could 
conceivably build up in more massive galaxies, but the availability of sufficient metal-free 
gas in these systems is unclear. Models for the evolution of the Pop III star formation rate 
density agree that some level of Pop III star formation persists down to z ~ 5 — 6 (see 
Figure 2), but disagree greatly on the average mass in young Pop III stars per galaxy. For 
example, Xu et al. (2016), Mebane et al. (2018), Liu & Bromm (2020b) and Skinner & Wise 
(2020) argue that close to reionization, the typical mass in young Pop III stars per galaxy 
is ~ 10? Mo, while Pallottini et al. (2014) and Sarmento et al. (2018) find much larger 
values of ~ 10? Mo. This difference stems from differences in how Pop III star formation in 
minihalos is treated in the models. If the process is strongly suppressed by LW radiation, 
baryonic streaming or the impact of limited numerical resolution, then there is less pre- 
enrichment accreted by more massive galaxies, making it more likely to find large clusters 
of Pop III stars in these systems. On the other hand, if Pop III star formation in minihalos 
is barely affected by early feedback, as in e.g. Skinner & Wise (2020), then there is greater 
pre-enrichment of the gas in massive galaxies, resulting in smaller Pop III cluster masses. 
Intriguingly, this suggests that a non-detection by JWST of any Pop III-dominated high 
redshift galaxies may allow us to constrain the efficiency of star formation in minihalos. 

Another way to reach the brightness necessary for JWST to detect small clusters of 
massive Pop III stars is gravitational lensing. For the example discussed above, we would 
need a magnification factor of around 100 in order to increase the brightness sufficiently. 
To detect individual Pop III stars at high redshift, even higher magnifications would be 
required. Rydberg et al. (2013) estimate that a magnification of around 3000 is required 
for JWST to detect a lensed 60 Mo star at z ~ 10. Such high magnifications are achievable 
close to lensing caustics (Windhorst et al. 2018; Diego 2019), but the disadvantage is that 
only a very small fraction of the high redshift Universe is magnified by such a large amount. 
The chance of actually detecting a highly lensed Pop III star therefore depends on the 
number of massive Pop III stars present in the high redshift Universe and hence on the Pop 
III SFRD. Unfortunately, even for values of this at the upper range of predictions from the 
literature (see Section 3.1.3), the probability of actually detecting a highly lensed Pop III 
star or star cluster is small (Rydberg et al. 2013; Diego 2019). Despite this, two possible 
candidate objects have recently been reported in the literature: a highly lensed star cluster 
at z — 6.629 (Vanzella et al. 2020) and Earendel, an extremely magnified single star at 
z — 6.2 (Welch et al. 2022), which Schauer et al. (2022) argue has a non-negligible chance 
of being a Pop III star. Whether either of these objects is truly primordial is currently 
unclear but will likely be resolved by JWST in the near future. 

Finally, if some fraction of Pop III stars are supermassive, with masses of M ~ 10° Mo 
and luminosities of L ~ 10? Lo, then these may be directly detectable by JWST with no or 
minimal lensing (Surace et al. 2018, 2019; Vikaeus et al. 2022). They may also be detectable 
by Euclid and the Nancy Grace Roman Space Telescope (RST), given magnifications of 
10 — 1000 (Vikaeus et al. 2022). However, we expect these stars to be rare objects (see 
discussion in Section 3.4), which strongly limits the number we expect to detect in a typical 
survey. For example, Vikaeus et al. (2022) show that with reasonable survey parameters 
RST will only detect ^ 10 supermassive stars, even if their number density at high redshift 
is close to the predicted upper limits, while Euclid will be unlikely to detect any. On the 
other hand, JWST should detect ~ 30 or more in the redshift range 7 < z « 10. 

It is also important to note that detecting Pop III stars in the early Universe is only part 
of the problem. Equally important is our ability to distinguish them from other high-redshift 
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objects. In the observer’s frame, Pop III stars will have very red colours owing to their high 
redshifts and the almost complete absorption of radiation shortwards of Lyman-a by the 
neutral IGM. This is true regardless of their rest-frame colours, but will be exacerbated 
in cases where pre-main sequence Pop III stars or highly-inflated massive stars contribute 
significantly to the observed flux (Mitani et al. 2019; Woods et al. 2021c). However, Pop II- 
dominated galaxies at the same redshift will also be very red, as will dusty starbursts at lower 
redshifts (Naidu et al. 2022). Pop III-dominated systems can potentially be distinguished 
from Pop II-dominated systems based on their broad-band colours, but doing this purely 
with the NIRCAM filter set requires an almost complete absence of nebular emission (i.e. 
a very high escape fraction of ionizing photons, fesc > 0.95; see Zackrisson et al. 2011). For 
galaxies with lower escape fractions, detection in the MIRI 560W and 770W filters is also 
necessary for distinguishing between Pop III and Pop II stars, limiting this method to only 
the brightest sources (Zackrisson et al. 2011; Trussler et al. 2022). 

A more promising route for distinguishing Pop III stars from other high redshift objects 
involves spectroscopy. Massive Pop III stars are unusually hot compared to metal-enriched 
stars (Schaerer 2002) and hence produce a much larger flux of photons with energies above 
54 eV that are capable of ionizing Het to Hei". It was therefore quickly realized that 
the recombination lines of Hell, particularly the 1640À line, could act as an important 
diagnostic of Pop III star formation (Oh et al. 2001; Tumlinson et al. 2001; Schaerer 2002). 
Unfortunately, detection of this line by itself does not prove that the source is a Pop III star, 
as Wolf-Rayet stars, X-ray binaries and black holes can also produce strong Hell emission 
(see e.g. Erb et al. 2010; Schaerer et al. 2019). A salutary example is provided by the z = 6.6 
Lyman-a emitter CR7 (Sobral et al. 2015). It shows a strong Hell 1640A emission line that 
was interpreted by some authors as being due to Pop III star formation (e.g. Pallottini et al. 
2015; Visbal et al. 2016). However, subsequent observations detected [OII] emission from 
this source (Bowler et al. 2017), which is inconsistent with this explanation, and it is now 
thought that CR7 is either a low luminosity AGN or a metal-poor starburst. A Pop III 
origin for Hell 1640À emission becomes more plausible if there is no associated emission 
from highly ionized metals, but even in this case a Pop II origin cannot be excluded (Katz 
et al. 2022). Altogether, the equivalent width of Hell recombination lines appears to be 
a promising diagnostic parameter, with Pop III stars producing larger widths than other 
astrophysical sources (Nakajima & Maiolino 2022), although the requirement to measure it 
accurately may again restrict the use to only the brightest sources. 


5.1.2. Pop IIl transients: supernovae and gamma-ray bursts. The high peak brightness 
of Pop III CCSN and PISN potentially make them much easier targets to observe than 
their massive star progenitors. The predicted peak luminosity of a Pop III CCSN can 
reach values as high as a few times 10!! Lo (Whalen et al. 2013a), around 10^ times larger 
than the luminosity of the single 1000 Mo star. Pop III PISNe are even brighter, with 
peak luminosities a factor of ten or more greater than for CCSNe (Whalen et al. 2013b). 
Although these peak luminosities persist for only brief periods, Pop III SNe remain brighter 
than their progenitors over a period of weeks to months. 

In a series of papers, Whalen and collaborators have examined the detectability of Pop 
III CCSNe, hypernovae, pulsational pair instability SNe and PISNe by JWST and other 
future facilities (Whalen et al. 2013a,b, 2014; Smidt et al. 2014, 2015). They show that Pop 
III PISNe will be bright enough to detect with JWST at redshifts 7 « z « 15 for hundreds 
of days. The brighter explosions may also be marginally detectable with the RST, but 
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will likely be too faint to detect with Euclid or in the Pan-STARRS or LSST sky surveys. 
Similar results are found for pulsational pair-instability SNe and hypernovae. CCSNe are 
harder to detect: Some explosion models are bright enough to be detectable with JWST at 
high redshift but others are not, and all are likely too faint to detect with the RST. 

Unfortunately, although many Pop III SNe will be bright enough for JWST to detect, 
we expect them to be rare because of their short duration, and so the expected number 
of detections is small. The Pop III SNe rate depends on the star formation rate and IMF 
of Pop III stars and hence is subject to the considerable uncertainties in both of these 
quantities that we have discussed previously. However, several recent studies broadly agree 
that we should expect around 3 — 7 x 1074 Pop III PISNe per JWST NIRCAM field of view 
per year (Hummel et al. 2012; Johnson et al. 2013; Magg et al. 2016; Hartwig et al. 2018a). 
The numbers for CCSNe are larger, but not by more than an order of magnitude (Magg 
et al. 2016). The chance of detecting any Pop III SNe in a dedicated survey is therefore 
tiny. Serendipitous detections remain plausible, although unambiguously identifying the 
detected objects as Pop III SNe without multi-epoch imaging may be problematic (Hartwig 
et al. 2018a). At redshifts close to the epoch of reionization, confusion with SNe produced 
by low metallicity Pop II stars will also be a major issue (Moriya et al. 2022). 

The other transients related to Pop III star formation that are detectable at cosmologi- 
cal distances are long-duration gamma-ray bursts (GRBs). The extremely high luminosities 
associated with these events makes them readily visible even at high redshifts, as demon- 
strated by GRB 090423, which has a spectroscopically-confirmed redshift of z ~ 8.2 (Tanvir 
et al. 2009) and GRB 090429B, which has a photometric redshift of z ~ 9.4 (Cucchiara et al. 
2011). Current gamma-ray telescopes such as Swift (Gehrels et al. 2004) are unlikely to 
detect GRBs at redshifts much higher than z ~ 9, because of limited sensitivity, and hence 
are mostly sensitive to GRBs produced by metal-enriched stars (Ghirlanda et al. 2015). 
However, future gamma-ray facilities such as THESEUS will expand the redshift range in 
which GRBs and their afterglows can be readily detected up to z ~ 20 (Amati et al. 2018), 
allowing us to probe redshifts at which Pop III sources will dominate. In view of this, it 
is clearly of interest to be able to make predictions for the number of Pop III GRBs that 
we expect these future facilities to detect. Unfortunately, this turns out to be even less 
well constrained than the number of Pop III SNe we expect at high redshift. As well as 
uncertainties due to the Pop HI SFR density and IMF, the number of GRBs that we expect 
to detect also depends on several other highly uncertain numbers (Bromm & Loeb 2006): 
the fraction of massive Pop III stars that are found in close binary systems, the fraction 
of these stars that are rapidly rotating, and the beaming factor (i.e. the probability that a 
randomly selected GRB happens to beam its radiation in our direction). Estimates of the 
number of Pop III GRBs that we should detect given a perfect telescope therefore range 
from values as high ~ 100 yr~! (Kinugawa et al. 2019; Lazar & Bromm 2022) to ones as 
small as 0.1 yr ^! (Toma et al. 2016). 


5.1.3. 21 cm global signal. The ground state of atomic hydrogen exhibits a hyperfine split- 
ting owing to the interaction between the proton and electron spins, with an energy differ- 
ence of 5.9 x 107° eV (corresponding to A ~ 21 cm) between the two states. The radiative 
transition between these states — the 21 cm line — is frequently used to determine the prop- 
erties of atomic hydrogen from the local Universe to very high redshifts. The global (i.e. 
sky-averaged) 21 cm signal is an important probe of the Universe at the redshifts when the 
first stars formed. When averaged on very large scales, the contribution of density fluctua- 
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tions cancel out and the strength of the line is determined solely by the spin temperature 
Ts of the transition, yielding an emission line if this is larger than the CMB temperature 
Tcmb, and an absorption line if it is smaller. The evolution of Ts is controlled by three 
physical processes: collisions with H, Ht and e , the absorption and emission of 21 cm 
photons, and the Wouthuysen-Field effect (i.e. hyperfine transitions brought about by the 
absorption and emission of Lyman-a photons; see Wouthuysen 1952, Field 1958). 

At redshifts immediately prior to the formation of the first stars, the interaction with the 
CMB dominates and Ts ~ Towns. However, once Pop III stars start to form, the Lyman-a@ 
photons that they produce allow the Wouthuysen-Field effect to become dominant, driving 
Ts towards the kinetic temperature of the gas. The timing of this transition therefore lets 
us learn something about when (and how quickly) Pop III stars started to form, and the 
information that it provides us on the temperature of the gas constrains the rate at which 
the high-redshift IGM is heated by X-ray sources (e.g. Fialkov et al. 2014; Pacucci et al. 
2014), its degree of ionization (which becomes more important at lower redshift; Park 
et al. 2020), and possibly also the coupling between dark matter and gas (e.g. Barkana 
2018; Muñoz et al. 2018; Liu et al. 2019). A measurement of the global 21 cm signal could 
also in principle test the existence of any radio background in addition to the CMB (e.g. 
Feng & Holder 2018; Ewall-Wice et al. 2018; Fialkov & Barkana 2019; Reis et al. 2020). 
Measurement of the high-z 21 cm line therefore provides stringent constraints on the early 
evolution of the cosmic star formation history and the properties of the first stars (Madau 
& Dickinson 2014), as reviewed in Furlanetto et al. (2006), Barkana & Loeb (2001), Fan 
et al. (2006a), Morales & Wyithe (2010), Pritchard & Loeb (2012), or Mesinger (2019). 

Experiments aiming at measuring the global signal at high redshifts include EDGES 
(Bowman et al. 2013), LEDA (Price et al. 2018), SARAS (Singh et al. 2018), PRIZM 
(Philip et al. 2019), or REACH (de Lera Acedo et al. 2022). The first tentative detection of 
this signal was reported by the EDGES collaboration (Bowman et al. 2018). Their signal is 
much stronger than expected and efforts to explain this have led to considerable theoretical 
speculation but few firm conclusions. In addition, we note that the true nature of this signal 
is highly debated (e.g. see Hills et al. 2018; Sims & Pober 2020). If it is of cosmological 
origin it represents one of our first direct constraints on the star formation process at a 
redshift of z ~ 17 as reflected by the impact of Lyman-a coupling and X-ray heating. 

The characteristics of early star formation can influence the global 21cm signal in a 
variety of additional ways. For example, the recovery time (i.e. the average duration for 
metal-enriched material, which has been expelled from the Pop III star-forming halo, to 
mix with the ambient material and recollapse to form Pop II stars; see Section 4.3.4) can 
influence the observable 21cm signal. Mirocha et al. (2018) and Magg et al. (2022b) find 
that short recovery times, resulting from weak supernova feedback, possibly due to Pop III 
star formation occurring in more massive halos or with low efficiency, result in a noticeably 
steeper global 21cm signal than models with long recovery times from strong feedback, 
characteristic of scenarios in which stars build up in low-mass halos or form in groups 
leading to mechanical energy and momentum input from multiple supernovae. The global 
21cm signal is also sensitive to the Pop III IMF. In a study focusing on the impact of 
Lyman-a radiation at high redshift and neglecting ionization and X-ray emission, Gessey- 
Jones et al. (2022) show that current and future observations can probe the characteristic 
stellar mass if the IMF is dominated by stars lighter than ~ 20 Mo. However, the expected 
signatures are relatively weak and comparable to the instrumental precision. If the Pop III 
IMF is dominated by very massive stars, we can no longer discriminate between different 
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IMF models. This may change once calculations including X-rays have been performed, 
which are thought to have a different dependence on the IMF compared to the Lyman-a 
emission. Detailed measurements of the 21 cm emission at high redshifts, once available, will 
open up many interesting pathways to study the birth of Pop III stars and the transition 
to Pop II star-formation. 


5.2. Low-redshift Universe 


The most stringent constraints on the properties and mass distribution of Pop III stars are 
likely to come from the study of extremely metal-poor (EMP) stars in the Local Group. 
These stars have a metallicity of [Fe/H] < —3 as expressed by the iron abundance relative 
to the solar value (for a review see Beers & Christlieb 2005). A common approach is to 
use large-area photometric or low-resolution spectroscopic surveys’ for a first estimate of 
the stellar metallicity (e.g. Ludwig et al. 2008; Caffau et al. 2011b). This information can 
then be used to put together a list of objects for high-resolution spectroscopic follow-up 
observations. Good targets are main-sequence turn-off stars, because they are numerous 
and bright, and they have chemical abundances that are essentially unaltered since the time 
of their formation (Pinsonneault 1997). The resulting stellar archeological surveys, such as 
those aiming at the Galactic halo and/or bulge (Christlieb et al. 2008; Schórck et al. 2009; 
Frebel 2010; Salvadori et al. 2010; Caffau et al. 2013; Starkenburg et al. 2017; Arentsen 
et al. 2020; Ishigaki et al. 2021; Li et al. 2022) or nearby satellite dwarf galaxies (e.g. Koch 
et al. 2013; Salvadori et al. 2015; Kirby et al. 2015; Roederer et al. 2016; Ji et al. 2016a,b, 
2023), can contribute to our knowledge of primordial stars in two key ways. 


5.2.1. Galactic archaeology and the high-mass end of the Pop III IMF. We can use the 
abundance pattern determined in EMP stars to infer the properties of the progenitor stars 
which provided the observed chemical enrichment. It is common practice to employ stellar 
evolution calculations combined with explosive nucleosynthesis models for this compari- 
son, both for individual stars (e.g. Ishigaki et al. 2014a; Bessell et al. 2015; Placco et al. 
2016; Aguado et al. 2018a, 2019) as well as for large samples (Fraser et al. 2017; Ishigaki 
et al. 2018). This is facilitated by libraries of modelled supernova yields (Heger & Woosley 
2010; Nomoto et al. 2006; Umeda & Nomoto 2008; Ishigaki et al. 2018), which typically 
depend on the stellar mass of the exploding star, the explosion energy and one or several 
parameters that quantify the mixing-and-fallback process. T'he latter cannot be simulated 
self-consistently in 1D supernova models and in principle require detailed 3D radiation- 
hydrodynamic calculations of the star-forming halo to determine (Magg et al. 2022c). An 
additional constraint on this process is provided by the fact that the supernova yields need 
to be diluted with metal-free gas in the surrounding of the exploding star in order to match 
the absolute metallicity of the observed star (Magg et al. 2019). 

'The observational data indicate that most EMP stars are rich in carbon in relative 
terms (Beers & Christlieb 2005; Frebel & Norris 2015). Of specific interest amongst these 
carbon-enhanced metal-poor (CEMP) stars are the so-called CEMP-no stars. Despite an 
excess of carbon relative to iron of more than [C/Fe] = 1 they show no enhancement in 


5For example, the Sloan Digital Sky Survey provides spectral data for several 10? stars (for DR12, 
see García Pérez et al. 2016), and Gaia for more than 10? sources (for DR3, see Gaia Collaboration 
et al. 2022). 
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neutron-capture elements, i.e. they have [Ba/Fe] < 0 (Frebel et al. 2005; Komiya et al. 
2007; Norris et al. 2013; Keller et al. 2014; Aguado et al. 2018b; Nordlander et al. 2019). 
The origin of the elemental abundance patterns in these stars is one of the key questions 
of early chemical evolution. Their dynamics appear to be inconsistent with them having 
gained their metals via mass transfer from a binary companion (Zepeda et al. 2022). Umeda 
& Nomoto (2002) propose that the abundance pattern of CEMP-no stars is the fingerprint 
of a so-called mixing-and-fallback or faint supernova (see also Ishigaki et al. 2014b; Chan 
et al. 2020). Such events have relatively low explosion energies, and consequently the inner 
layers of the star, containing in particular iron, are suggested to fall back on the compact 
remnant, such that only the outer layers, containing carbon and other light elements, are 
ejected and become mixed with ambient material. 

Assuming that the oldest and most metal-poor stars in the Galaxy have been supplied 
with heavy elements by only one or at most two supernova explosions — a common assump- 
tion but one that may not be true for all EMP stars (Hartwig et al. 2018b, 2019; Welsh et al. 
2021) — one finds that the relative abundances of heavy elements are most consistent with 
core collapse supernovae from Pop III stars in the mass range 20 — 40 Mo (e.g. Bonifacio 
et al. 2015; Caffau et al. 2011a, 2012; Cooke & Madau 2014; Frebel et al. 2005; Iwamoto 
et al. 2005; Joggerst et al. 2009, 2010; Norris et al. 2013; Keller et al. 2014; Lai et al. 2008), 
although a few examples exist of stars that are more consistent with enrichment by hyper- 
novae (Placco et al. 2021; Skúladóttir et al. 2021) or PISNe (Aoki et al. 2014; Salvadori 
et al. 2019). The conclusion that most Pop III SNe were CCSNe is consistent with the work 
of Mapelli et al. (2006), who compute the expected density of intermediate-mass black holes 
in the Galaxy by assuming they are the relics of higher-mass Pop III stars. This allows 
them to derive an upper limit from the non-detection of ultra-Iuminous X-ray sources in the 
Milky Way. The approach has been extended by de Bennassuti et al. (2014, 2017) who also 
include the Galactic metallicity function and the abundances of EMP stars. Again these 
data appear inconsistent with pair-instability supernovae being important (although see 
Salvadori et al. 2019 for a contrary view). Taken together, these studies suggest that most 
Pop III stars forming through the standard pathway outlined in Section 3.2 are not likely 
to reach masses significantly above ~ 40 Mọ, with more massive objects being extremely 
rare and possibly being the result of collisions in very dense clusters (see Section 5.3). 

Finally, we note that complementary information on the metal yields produced by the 
first generation of stars is potentially provided by observations of extremely metal poor 
foreground absorption systems in the spectra of distant quasars. The most metal-poor of 
these systems have metallicities comparable to EMP stars (e.g. Cooke et al. 2017; Robert 
et al. 2022; Welsh et al. 2022). If we assume that each of these systems has been enriched 
by either a single Pop III supernova or a small number of such SNe, then we can use the 
observed metal abundances to infer their properties. In general, the results we obtain from 
this procedure are consistent with those coming from EMP stars. Most of the systems 
show elemental abundances consistent with enrichment primarily by CCSNe (Welsh et al. 
2019), although a recent example has been found of a cloud in the broad-line region of 
a z = 7.6 quasar that has an abundance pattern that appears to indicate a PISN origin 
(Yoshii et al. 2022). This supports the idea that most Pop III supernovae are CCSNe, with 
PISNe occurring only rarely. 


5.2.2. Potential for direct detection of Pop III stars and low-mass end of the IMF. Stellar 
evolution models (e.g. Kippenhahn et al. 2012; Maeder & Meynet 2012) tell us that low-mass 
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Figure 12 


Properties of the Milky Way and its satellite galaxies. The horizontal line in each column 
indicates the stellar mass as inferred from observations. The filled circles indicate the probability 
that the systems contains at least one low-mass Pop III survivor, with the color encoding the 
median fraction. The figure is based on data from Magg et al. (2018). 


Pop III stars with M < 0.8Mo6 must have survived until the present day. Furthermore, 
models of the surface pollution of Pop III stars by metals from the ISM suggest that this 
is not an important effect and that Pop III stars around at the present day would still be 
spectroscopically identifiable as such (Frebel et al. 2009; Tanaka et al. 2017). Therefore, 
if these stars exist, there is a chance to directly detect some of them in current or future 
stellar archeological surveys. One can use semi-analytic models to predict the most likely 
locations to find these Pop III survivors. Whereas the predicted numbers of survivors 
are highest in the Galaxy, specifically in the bulge, they are very difficult to find there 
because of the presence of so many higher metallicity stars. It is therefore better to search 
in environments that did not have continuous star formation over secular timescales, but 
instead only experienced an initial burst of star formation at high redshift and then remained 
quenched during subsequent cosmic evolution. Consequently, the lowest-mass satellites of 
the Milky Way are expected to contain the largest fraction of Pop III survivors compared to 
their Pop II stellar content, and so they are the most promising targets for future attempts to 
find Pop III survivors with high-resolution, high signal-to-noise spectroscopic observations 
(Magg et al. 2018). Figure 12 shows predictions of the expected fraction of genuine low-mass 
Pop III stars in Local Group galaxies together with their total stellar mass. We see that 
the lowest-mass satellites are the best places to look for these stars, although it is possible 
that some of these systems will be too small to contain any Pop III survivors. 

Even non-detections allow us to put stringent limits on the low-mass end of the Pop III 
IMF. For example, Hartwig et al. (2015) estimate the expected numbers of low-mass Pop III 
stars in the Galactic halo based on semi-analytic models of the early star formation history 
in Milky Way-like halos. They conclude that if no such object is found in surveys with 
sample sizes of 4 million stars then we can exclude the existence of low-mass Pop III stars 
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with masses below 0.8 Mo with high significance. On similar grounds, Salvadori et al. (2007, 
2010) develop a detailed model of the metallicity distribution function of metal-poor stars in 
the Galactic halo and suggest that Pop III stars should be more massive than 0.9 Mo when 
matching their predictions with the observational data. A different conclusion was reached 
by Tumlinson (2006, 2010), who argues that the current abundance measurements are not 
precise enough to distinguish between different models of the primordial IMF based on 
studying the chemical evolution during the early build-up of the Galaxy. However, he also 
suggests that characteristic masses of a few 10 Mo provide a better fit to the available data 
than masses of 100 Mo or above, consistent with the supernova yields mentioned above. 
Finally, Rossi et al. (2021) argue that our failure to find any Pop III stars in nearby ultra- 
faint dwarf galaxies already allows us to exclude Pop III stars with M < 0.8Mo. Altogether, 
the prospects of finding surviving low mass first stars in our immediate neighborhood are 
highly exciting, but at present no firm conclusions about the low-mass end of the Pop III 
IMF can be drawn. This may change with the availability of new generations of multi-object 
spectrographs, such as 4MOST or WEAVE (e.g. Dalton et al. 2012; Feltzing 2016; Feltzing 
et al. 2018), which will allow for galactic archeological surveys of unprecedented size and 
depth that can help to resolve this issue. 


5.3. Gravitational waves 


The detection of gravitational waves with extremely precise interferometers (advanced 
LIGO, Aasi et al. 2015, advanced Virgo, Acernese et al. 2015, or KAGRA, Aso et al. 
2013, all of which are currently operational, and the upcoming LISA, Amaro-Seoane et al. 
2023, DECIGO, Kawamura et al. 2011, and TianQin Luo et al. 2016 experiments in space 
or the Einstein Telescope on the ground, Punturo et al. 2010) opens up a new observational 
window to the Universe. The signals detected so far come from the gravitational inspiralling 
and merging of two stellar-mass black holes, two neutron stars, or pairs of a black hole and 
a neutron star. Whereas LIGO/Virgo/KAGRA can detect the mergers of compact objects 
with masses up to several hundred solar masses out to z ~ 1 (e.g. Abbott et al. 2019), the 
planned observatories (such as LISA, TianQin, and DECIGO, or the Einstein Telescope) 
have the potential to detect considerably more massive objects out to higher redshifts. 
There are two primary ways to infer information about primordial stars from these data. 
On one hand, we can look at the properties of individual measurements, which gives us in- 
formation about the total mass of the system, the mass ratio of the two constituents, and 
their spin. This provides constraints on the mass spectrum and multiplicity of Pop III stars 
and on possible formation and evolution pathways. On the other hand, we can study the 
rate of these events (typically normalized to a cosmic volume of 1 Gpc?) as function of red- 
shift, and compare the results with our understanding of the cosmic star-formation history. 
However, the astrophysical interpretation of these data is highly challenging for several rea- 
sons. First, the detected gravitational wave signal from the merger of two compact objects 
can occur many Gyr after the formation of the original stellar binary system. To link these 
two events requires complex binary evolution calculations that reproduce the correct stellar 
structure in the MS and post-MS phases and that reliably model mass transfer and the 
change of orbital parameters over secular timescales (see e.g. Belczynski et al. 2002, 2016; 
Eldridge & Stanway 2016; de Mink & Mandel 2016; Mandel & de Mink 2016; Marchant 
et al. 2016, 2021; Stevenson et al. 2017; Mapelli et al. 2017; Spera et al. 2019). Second, to 
explain the current data, the mergers of stars in a dense cluster environment may be neces- 
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sary. The corresponding stellar dynamical processes contribute additional stochasticity and 
uncertainty to the problem (e.g. Portegies Zwart & McMillan 2002; Portegies Zwart et al. 
2022). And finally, the signal coming from genuine Pop III stars is much weaker than the 
rates expected from metal-enriched binaries, simply because the total number of metal-free 
stars is many orders of magnitude smaller, which can be see for example by comparing 
the corresponding cosmic star formation rate densities in Figure 2, and which cannot be 
compensated by differences in the IMF (Figure 11). 

The combined data of the first three observing runs of LIGO, Virgo, and KAGRA 
(Abbott et al. 2021a,b) demonstrate that the mass spectrum inferred for the primary and 
secondary merger component is very wide and contains structure beyond a simple power 
law with a sharp high-mass cut-off. The mass distribution of the primary shows a first peak 
at ^ 10Mg and a secondary one at 30 — 35 Mo. Although the majority of the systems 
detected have primary masses smaller than 45 Mo, the mass distribution clearly extends 
beyond ~ 65 Mo and reaches into the so-called pair-instability gap, at least for metal-rich 
stars (e.g. Farmer et al. 2019; Woosley & Heger 2021; see also Section 4.3.4). Pair-instability 
and pulsational pair-instability supernovae lead to the complete disruption of the star and 
leave no remnant behind. Because metal-rich stars experience considerable mass loss due to 
strong winds or pulsational instabilities prior to the explosion, we expect that no black hole 
in the mass range between about 50 Mo to 130 Mo should form, whereas for zero-metallicity 
stars, which experience no mass loss, the pair-instability gap covers masses from ~ 90 Mo 
to — 260 Mç (however, with considerable uncertainties, see e.g. Heger & Woosley 2002; 
Woosley & Heger 2015; Farmer et al. 2020; Costa et al. 2021; Vink et al. 2021). 

As we expect that most of the observed gravitational wave detections come from the 
remnants of metal-enriched rather than zero-metallicity stars (e.g. Hartwig et al. 2016; 
Belczynski et al. 2017; Tanikawa et al. 2022), these findings cast doubts on the widespread 
occurrence of pair-instability supernovae, and it has been argued that stellar collisions in 
a dense cluster environment provide an alternative pathway to produce the high-mass end 
of the observed black hole mass spectrum (e.g. Portegies Zwart & McMillan 2000, 2002; 
Portegies Zwart et al. 2004; Giersz et al. 2015; Mapelli 2016; Di Carlo et al. 2019, 2021; Liu 
& Bromm 2020a; Dall' Amico et al. 2021; Wang et al. 2022; Costa et al. 2022). At the low 
mass end of the distribution, the data indicate a lack of low-mass black holes below ~ 6 Mo 
with high significance. This is larger than the masses reported for black hole candidates 
in Galactic binaries (reaching down to ~ 3 Mọ, see Thompson et al. 2019) and potentially 
hints towards a more top-heavy IMF of the progenitor stars. Most systems have mass ratios 
q close to unity, but there are some for which q « 1 is inferred with high significance. 

Another feature of these reported gravitational wave detections is that about one third 
of the systems exhibit spins which are possibly inclined with respect to the orbital angular 
momentum and which show signs of precession (see also Callister et al. 2022). This provides 
additional support for a collisional origin of these systems. The inferred merger rates for 
binary black holes have a 90% credibility interval of 18 — 44 Gpc ? yr~'; those for binary 
neutron stars lie in the range of 10 to 1700Gpc ? yr ! (Abbott et al. 2021b). These 
values are still highly uncertain, but the precision will increase with time as the number 
of detected gravitational wave events continues to go up. The data furthermore indicate a 
merger rate that grows with redshift as (1+ z)^, with « only marginally constrained in the 
range from 1.1 to 4.6, at least for z < 1. Theoretical predictions of gravitational rate event 
rates, specifically from the remnants of Pop III stars (Kinugawa et al. 2014; Hartwig et al. 
2016; Schneider et al. 2017; Belczynski et al. 2017, 2022; Dayal et al. 2019; Neijssel et al. 
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2019; Tang et al. 2020; Ng et al. 2021; Tanikawa et al. 2021, 2022), are roughly consistent 
with these numbers, but more stringent constraints clearly require larger statistics (for a 
comprehensive overview of this aspect, see Mandel & Broekgaarden 2022). 


6. Summary 


Summary points 


1. The formation of the first stars begins at redshifts of z ~ 30 and is likely to continue 
to z ~ 5, dominating the cosmic star formation rate density until z ~ 15, at which 
point the build-up of metal-enriched Pop II stars takes over (Section 3.1). 

2. Stellar birth in the high-redshift Universe shares the same complexity as present- 
day star formation. Primordial gas is highly susceptible to fragmentation and we 
expect the first stars to form in binary or higher-order multiple stellar systems with 
a wide range of separations, a flat mass ratio distribution, and a roughly thermal 
spread of orbital eccentricities (Section 3.2.2). 

3. Most primordial stars build up in isolated high-redshift halos, which have not been 
affected by stellar feedback. We call this the standard pathway of Pop III star 
formation. We expect a wide distribution of masses, potentially ranging from the 
sub-solar regime up to several hundred solar masses. The initial mass function 
(IMF) is approximately logarithmically flat, which makes it top-heavy compared to 
the Milky Way values (Section 3.2.3). 

4. Feedback from the first massive Pop III stars plays a major role in regulating subse- 
quent star formation. Photoionization and supernovae dominate within individual 
halos, whereas the Lyman-Werner and X-ray backgrounds built up by Pop III stars 
and their remnants govern the evolution on cosmic scales (Sections 4.2 and 4.3). 

5. Pop III supernovae produce metals that enrich gas associated with their host halo 
(“internal enrichment") but can also reach neighboring halos (“external enrich- 
ment"). The time it takes Pop II stars to form from metal-enriched gas remains 
poorly constrained and depends on the halo mass, the supernova energy, and de- 
tails of the local environment. It is plausible that some fraction of the observed 
extremely metal-poor stars in the Milky Way and its satellites were enriched by 
single supernovae (Section 4.4). 

6. If they exist, Pop III stars with masses below ~ 0.8 Mọ will have survived until 
today and should be detectable in surveys of extremely metal-poor stars. However, 
none has yet been found. At the high-mass end, observations point towards Pop 
III stars that end their lives in core-collapse supernovae, implying a mass range of 
20 — 40 Mo. Evidence for more massive stars is sparse (Sections 5.1 and 5.2). 

7. Under very rare and extreme conditions, the formation of supermassive stars is 
possible. Subjected to accretion rates of M ~ 0.1 Mo yr! or higher, they can 
grow to masses of several 10° Mo, before the general-relativistic instability kicks in 
and triggers the collapse to a black hole. These objects can be the seeds for the 
supermassive black holes observed in high-redshift quasars. Because the infalling 
gas is expected to fragment, run-away collisions in dense and deeply embedded 
stellar system are thought to be important in this process (Section 3.4). 
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Future Issues 


il, 


From a theoretical point of view, better understanding the impact of stellar feedback 
is a key challenge of current research into Pop III star formation. Specifically, the 
question of how photoionizing radiation escapes from the immediate vicinity of the 
star remains uncertain. On large scales, the relative importance of Lyman-Werner 
feedback, X-rays and streaming velocities remains unclear, with several recent stud- 
ies coming to highly different conclusions. Overall, it is likely that feedback has a 
negative effect on the Pop III star formation rate, but its influence on the mass 
distribution and the multiplicity of Pop III stars remains undetermined. 

On the observational side, no genuine low-mass Pop III star has been found so far 
in stellar archeological surveys. With the increasing sample size provided by future 
surveys, it will become possible to either confirm or rule out the existence of Pop 
III stars with masses below ~ 0.8 Mo with high statistical significance. 

Halos forming Pop III stars are unlikely to be observable with JWST without the 
aid of gravitational lensing, unless they have higher masses and star formation 
efficiencies than current models predict. Supermassive Pop III stars or extremely 
strongly lensed systems should be bright enough to see, but will be very rare, and so 
it is plausible that none will be detected (Section 5.1.1). Pop III supernovae should 
be visible with JWST, but are also very rare. Nevertheless, serendipitous discoveries 
remain a realistic possibility. Pop III gamma-ray bursts (GRB) are unlikely to be 
detected using current facilities, but will become accessible by next-generation GRB 
telescopes (Section 5.1.2). Much work still needs to be done to understand how to 
relate detections (or non-detections) of high-redshift Pop III stars to constraints on 
models of Pop III star formation. 

Detection and characterization of the global 21 cm signal at high redshift will allow 
us to put stringent constraints on the properties of high-mass Pop III stars as well 
as the early evolution of the cosmic star formation rate density. 

'The detection of gravitational waves from merging binary black holes in the mass 
range forbidden by pair-instability models provides additional support for the im- 
portance of stellar collisions, possibly occurring in the star clusters that are a natural 
consequence of the standard Pop III formation pathway (Section 5.3). Our ability 
to use these events to constrain the physics of Pop III star formation will improve 
as our sample size grows. 
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Appendix: Average stellar properties and production rates of ionizing and 
non-ionizing photons 


The stellar properties discussed in Section 4.1.2 and summarized in Figure 8 can be used 
to estimate the number of ionizing and non-ionizing photons emitted by stars of different 
mass and metallicity. We focus our discussion on the standard pathway of Pop III formation 
(Section 3.2), that is on stars which have experienced small to moderate accretion rates, 
and which therefore enter the main sequence (MS) compact and hot. Our approach is based 
on evolutionary tracks of individual stars that are computed by self-consistently solving the 
equations of stellar structure, the production and transport of energy in the interior of the 
star, and the loss of mass and radiation through its surface as a function of time in spherical 
symmetry (for further details, see the book of Kippenhahn et al. 2012). We consider both 
non-rotating stars as well as stars that spin with 40% of the break-up velocity (e.g. Maeder 
& Meynet 2012). The data presented here are obtained with the Geneva code (Eggenberger 
et al. 2008), and they are derived from stellar evolution models provided by Ekström et al. 
(2012), Georgy et al. (2013), Groh et al. (2019), Murphy et al. (2021a,b), and Martinet et 
al. (in prep.). We note that other approaches lead to very similar results. Specifically, we 
have compared with stellar evolution models obtained with the SEVN code (Spera et al. 
2022) and found differences of at most a few percent (M. Mapelli & T. Hartwig, private 
communication). 
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Figure 13 


Stellar parameters as function of time for selected non-rotating stars of metallicity Z = 0 (blue) 
and Z = Ze (red) with masses of 2Mo, 5Mo, 20 Mo and 120 Mọ (different line strength). From 
left to right: stellar mass, bolometric luminosity, effective surface temperature, and stellar radius. 
The tracks have been computed with the Geneva stellar evolution code by Murphy et al. (2021a,b) 
and Ekstróm et al. (2012). 


It is clear that key stellar parameters, such as mass, bolometric luminosity, effective 
temperature and radius, change during the main-sequence (MS) and post-MS evolutionary 
phases. We illustrate this in Figure 13, which shows these parameters as function of time 
for non-rotating stars of mass 2Mo, 5 Mo, 20 Mç and 120 Mç and zero as well as solar 
metallicity. When looking at the left-most panel, we notice that none of the zero metallicity 
Pop III stars experiences mass loss, whereas high-mass solar metallicity Pop I stars shed 
significant amounts of material in the post-MS phase due to strong line-driven winds and/or 
pulsational instability. The bolometric luminosity remains roughly constant during the 
main sequence, but can easily rise by a factor of ten or more in the giant phase, which 
is most noticeable in the second panel for stars of intermediate masses (M = 2Mo and 
M =5Mo). Independent of metallicity the effective surface temperature drops sharply 
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(third panel), while the stellar radius increases enormously (right-most panel) in the giant 
phase of post-MS evolution. 

We use these stellar evolution calculations to compute the spectral energy distribution 
at each point in time. For simplicity, we assume pure blackbody radiation and neglect in- 
dividual line emission and absorption features (see Schaerer 2002, for a critical assessment 
of this approach). We then integrate or take the average of the corresponding radiative 
flux over the entire MS and post-MS evolutionary phase. This allows us to ask how many 
ionizing and non-ionizing photons Pop III stars of different masses produce over their life- 
time. Specifically, we are interested in photons in the Lyman-Werner bands (with energies 
in the range 11.2eV € hv < 13.6eV), as these are important for the dissociation of H> 
(Sections 4.2.2 and 4.3.2), and in photons that can ionize hydrogen (H: hv > 13.6eV), 
neutral helium (Hel: hv > 24.6eV), and singly ionized helium (Hell: hv > 54.4eV). The 
results are presented in Figure 14, which not only considers the zero and solar metallicity 
case, but also values in between, roughly corresponding to what we expect for the Large and 
Small Magellanic Clouds (Z = 0.42 Zo and Z = 0.14 Zo, respectively) and for low-mass 
low-metallicity dwarf galaxies, such as I Zwicky 18 (Z = 0.03 Zo). 

We consider models without and with stellar rotation at the 40% break-up level. The 
non-rotating models tend to give slightly larger fluxes at all frequencies, but this is only 
noticeable in the low-mass regime. At the high-mass end, the difference between non- 
rotating and rotating stellar evolution models is small (see also Maeder & Meynet 2012). 
In the top row we provide the photon flux averaged over the stellar lifetime, in the middle 
row we present the total number of photons produced by each star during the MS and 
post-MS evolution, and in the bottom row the corresponding total number of photons per 
stellar baryon. The data entering Figure 14 are made available for further use in Table 
A1:, which we provide for all metallicities and both levels of stellar rotation in electronic 
form.” Here we only reproduce the entries for zero metallicity and non-rotating stellar 
models. Note that the lifetime-averaged data provided in the Table are close to the zero- 
age main sequence (ZAMS) values. This has two reasons. First, the changes of mass, 
luminosity, effective surface temperature, and radius are small as the star evolves through 
the MS. Second, the time the star spends on the post-MS (during which these parameter 
can vary by orders of magnitude, as illustrated in Figure 13) is comparatively short, and 
so this evolutionary phase does not contribute much to the overall photon budget. This is 
especially true for the high frequency range relevant for ionization, as the effective stellar 
temperature drops significantly on the post-MS and with it the corresponding UV flux in 
the blackbody approximation adopted here. 

We see that for high-mass stars above ~ 10 Mo the average fluxes in the LW band and 
for H ionization at different metallicities are comparable to each other within a factor of 
order unity. This may seem surprising at first glance, given the differences in effective surface 
temperatures (see Figures 8 and 13), but can be explained by subtle differences in stellar 
lifetime and post-MS evolution largely off-setting each other. When considering the hard 
UV emission necessary for helium ionization, the impact of the much higher effective surface 
temperatures of Pop III stars becomes clearly noticeable. These stars produce an order of 
magnitude more photons than their low-metallicity or solar metallicity counterparts, thus 


SThe full dataset is available from ARAA or from the authors (at heibox.uni- 
heidelberg.de/f/6b5b3fcbd3974fb98d50/) as a space-separated ASCII table, which can be read with 
the python astropy.io.ascii.read command. 
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Figure 14 


Stellar emission properties for five different metallicities (dark blue: Z = 0, light blue: 

Z = 0.03 Zo, green: Z = 0.14 Zo, orange: Z = 0.42 Zo, red: solar metallicity) and masses 
covering the range 0.8 Mç € M < 120 Mç for Z > 0.14 Zo, 1.7 Mç € M < 120 Mọ for 

Z < 0.03 Zo, plus four additional models 180 Mç < M < 500 Mọ for Z = 0. We focus on the 
Lyman-Werner band (11.2eV < hv < 13.6 eV) and on photons that can ionize H (hv > 13.6 eV), 
He (hv > 24.6 eV), and Het (hv > 54.4eV). The top row depicts the number of photons emitted 
per second averaged over the entire MS and post-MS lifetime. The middle row provides the total 
number of photons produced, and the bottom row gives the total number of photons per baryon in 
the star. The figure is constructed from data obtained from Ekstróm et al. (2012), Georgy et al. 
(2013), Groh et al. (2019), Murphy et al. (2021a,b) and Martinet et al. (in prep.). The tabulated 
values are available as machine readable table from ARAA or directly from the authors at 
heibox.uni-heidelberg.de/f/6b5b3fcbd3974fb598450/. 


opening up a potential observational window to discriminate between the contributions of 
Pop III and Pop II stars to early reionization by looking for differences in the volume- 
averaged rate of H and He ionization in the early Universe (Section 5.1.3). 
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The total number of photons per stellar baryon (bottom row of Figure 14) is a highly 
useful quantity when estimating the impact of star formation on the galactic and inter- 
galactic environment in numerical simulations and theoretical models. Combined with 
information about the stellar initial mass function (IMF) and linked to the cosmic star- 
formation rate density (see Figure 2) it gives an indication of the volume-averaged photon 
flux as function of redshift, which in turn is a prerequisite for estimating the evolution of 
the cosmic radiation background in different frequency bins. We compute the numbers for 
a representative population of zero-metallicy stars and compare to solar-metallicity stars 
in the Milky Way today. For Pop III we assume a logarithmically flat mass spectrum in 
the range 0.08 < M/Mo < 500 (see Section 3.2.3), whereas we adopt a (Kroupa 2002) 
multi-component power-law IMF with 0.08 < M/Mo x 120 for the solar neighborhood (see 
Section 4.4). The resulting population averaged photon numbers per baryon are listed in 
Table A2: in log units. We find that the more top-heavy IMF of Pop III stars and their 
larger mass range makes the difference in number of ionizing photons emitted per baryon 
for the population average even larger than was already noticeable when comparing indi- 
vidual stars. Altogether, we expect over 52000 hydrogen ionizing photons per baryon for 
the typical stellar population in the early Universe, and the numbers are about 19000 and 
1100 for the ionization of He and He*. In contrast, at the present day we have on average 
fewer than 500 photons per stellar baryon that can ionize H, 30 for He, and none for He”. 
This enormous discrepancy explains why primordial stars are so efficient in ionizing their 
environment and why they are key drivers of early cosmic reionization. 


Table A2: Population averaged number of photons per baryon (in log units) 


population IMF (Nb (Nb) (Nb) (Nó (NB a) 
Pop II log-flat 5.31 4.11 4.72 4.29 3.05 
Pop I Kroupa 6.11 2.67 2.69 1.50 0 


Column description: population = metal-free Pop III or solar-metallicity Pop I stars M IMF = stellar 
initial mass function, either logarithmically flat or (Kroupa 2002) multi-component power-law model M 
(N°) = decadic logarithm of the population averaged total number of photons emitted per stellar baryon 
= (New) = decadic logarithm of the population averaged number of LW photons emitted per stellar 
baryon M (N) = decadic logarithm of the population averaged number of H ionizing photons emitted 
per stellar baryon M (Neer) = decadic logarithm of the population averaged number of He ionizing 
photons emitted per stellar baryon M (Nien) = decadic logarithm of the population averaged number of 


Het ionizing photons emitted per stellar baryon 
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